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I. INTRODUCTION 



Every finite-impulse response (FIR) filter has two distinct properties [Ref. 5] ; 
first, it is always stable; second, if it is not causal, it can always be made to be causal 
by introducing a finite delay. FIR filters can be designed so that their frequency 
responses have an exactly linear phase characteristic. FIR filters with linear phase are 
important in applications like speech processing and data transmission, because in 
these applications a nonlinear phase filter is harmful. The symmetry property of the 
linear phase FIR filters, however, helps reduce the number of coefficients by nearly one 
half resulting in substantial computational savings. 

The various filter realizations, or structures that are frequently considered are the 
direct form, the cascade form, the parallel form, and the lattice form. The lattice form 
realization is of particular interest because of its superior numerical performance, and 
modularity in the structure. The operation of a lattice filter is completely described by 
specifying the sequence of reflection coefficients that characterize the individual stages 
of the filter. 

Conventionally an adaptive filter is composed of a tapped delay line or 
transversal structure with adjustable coefficients or weights and an adaptive algorithm 
which updates the coefficients continuously based on some performance criterion. The 
design of a fixed coefficient filter is based on the prior knowledge of both signal and 
noise. Adaptive filters, on the other hand, have the ability to adjust their own 
parameters automatically, and their design requires little, or no a priori knowledge of 
signal or noise characteristics [Ref. 14]. However, the designer has to choose the order 
of the filter and the type of the algorithm. Also the adaptive filter usually requires a 
large initial transient time (i.e., the initial filter convergence period). 

The least-mean-square (LMS) adaptive algorithm minimizes the mean square 
error c(k) by recursively altering the filter coefficient vector 3(k) at each sampling 
instant. The original Widrow-Hoff LMS algorithm is 3(k+ l) = 3(k) + 2pe^k)X(k), 
where ft is a convergence factor controlling stability and rate of adaptation [Ref. 15]. 
The algorithm is based on the method of steepest descent, moving b(k) in proportion 
to the instantaneous gradient estimate of the mean square error. The successive 
orthogonalization provided by the lattice offers advantages in adaptive convergence 
rate which cannot be achieved with tapped delay lines [Ref. 17,18]. 
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Recently, Prony's method has been applied to the estimation of spectral lines in 
noise [Ref. 11]. It has been shown that the Prony's method yields better spectral 
estimates than a companion spectral estimation technique, called Pisarenko's method 
[Ref. 10]. Also it has been established that the filter structure involved in Prony's 
method has a iinear phase property which is not the case for Pisarenko's method [Ref. 
11]. 

The thesis investigates the application of lattice structures in Prony's method of 
spectral line estimation. The complete solution to Prony's method consists of three 
steps: (i) representing a given process of M sinusoids in noise in terms of complex 
exponentials, (ii) finding roots of a symmetric polynomial, and (iii) estimating the 
frequency, phase and amplitude information. In this thesis, however, we have 
emphasized only the frequency estimation problem. Here we derive an adaptive lattice 
structure to realize linear phase transfer functions which will be used to estimate the 
sinusoidal frequencies as part of Prony's method. The scope of the thesis consists of 
obtaining a linear phase lattice structure, developing a least mean square (LMS) error 
based adaptive algorithm, and testing the lattice structure and the algorithm by means 
of computer simulation. 

The thesis is organized as follows. In Chapter II, we present a brief review of 
Prony's method of representing a given process in terms of a set of complex 
exponentials, and then address the problem of estimating spectral lines using this 
method. We show that the original Prony's method has to be modified slightly in 
order to apply it to the spectral line estimation problems. In Chapter III, we discuss 
the basic concepts of the linear phase FIR filter and the related lattice structure. We 
show three examples of obtaining the lattice structures for both linear phase and non- 
linear phase FIR transfer functions (Appendix A). In Chapter IV, we briefly discuss 
the least-mean-square (LMS) adaptive algorithm that results from attempting to 
minimize the mean square error and present a summary of some LMS algorithms for 
lattice that have been reported in the literature. Also included are some computer 
simulation results. In Chapter V, we present the derivation of a new LMS based FIR 
adaptive lattice algorithm and extend it to the linear phase case as well. The new 
algorithm is then used in the estimation of spectral lines in white noise. Results of 
simulation are included. 
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II. PRONY'S SPECTRAL LINE ESTIMATION 



A. INTRODUCTION 

In its original form, Prony's method analyzes processes involving simple, or 
damped harmonics using complex exponential functions as coordinate functions. On 
the other hand, the well known Fourier analysis consists of representing a given 
process in terms of a set of sine and cosine functions. Recently, Prony's method has 
been applied to the estimation of spectral lines in noise [Ref. 11], It has been shown 
that the Prony's method yields better spectral estimates than a companion spectral 
estimation technique, called Pisarenko's method, in terms of bias, spurious responses, 
and the computational complexity [Ref. 10]. 

In this chapter, we present a brief review of the Prony's method of representing a 
given process in terms of a set of complex exponentials, and then address the problem 
of estimating spectral lines using this method. We show that original Prony's method 
has to be modified slightly in order to apply it to the spectral line estimation problems. 

B. PRONY'S METHOD 

Consider that the given process, x(k), consists of n distinct sinusoids, then x(k) 
can be approximated by an expression of the form 



n 

x(k) =. J£ [A. costOjk + B ; sinco ; k] 



( 2 . 1 ) 



where the co/s are sinusoidal frequencies. The above approximation can be considered 
as a special case of an exponential approximation given by [Ref. 29] 



m 

V 



x(k) = Z C 

i= 1 1 






( 2 . 2 ) 



where m = 2n and the a/s are identified as C 0 j and -tOj. The values of co ; can be 
estimated, by Prony's method, assuming that the data are known at k = 0, 1, ..., N-l. 
From Eq.(2.2), we can obtain the following set of equalities 
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C 1 + C 2 + "• + C m = x (°) 

Cj e* a l + c 2 ei a 2 + ... + c m ei a m = x(l) 

C x eJ' 2a l + C 2 eJ 2a 2 + ... + C m ei 2a m = x(2) (2.3) 



Cj e^ N * l ) a l + C 2 e^ N '^ a 2 + ... + C m e^'^m = x(N-l) 

Now we have a set of N equations with 2m unknowns, namely, C. and a ; (i= 1, 2, ..., 
m) which can be solved exactly if N = 2m, or approximately by the method of least 
squares if \>2m. Also note that the N equations are nonlinear in the exponential 
terms e^ a i. Let e^ a i, i= 1, 2, ..., m, be the roots of the equation [Ref. 29]. 



e jma + a ^ e j(m-l)a + a ^ e j(m-2)a + _ + a m .i eja + a m = 0 (2.4) 

In order to determine the coefficients a ; (i = 1,2, ..., m), we multiply the first equation 
in Eq.(2.3) by a m , the second equation by ..., the equation by ctj and the 

(m+ 1)^ equation by 1, and add the results. In this way, we can obtain the (N-m) 
linear equations. From Eq.(2.2), we can obtain the following set of equalities 

x(m) + x(m-l) a x + x(m-2) a 2 + ... + x(0)a m = 0 
x(m+l) + x(m)a l + x(m-l) « 2 + ... + x(l)a m = 0 



(2.5) 



x(N-l) + x(N-2)a i + x(N-3)a 2 + ... +x(N-rn-I)a m = 0 
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Since the data x(k), k = 0, 1, 2, N-l, are known, this set of equations can be solved 
for the m a's if N ^ 2m. After the a's are determined, e ]a i, i = 1, 2, ..., m, are found 
from Eq.(2.4). The set of equations in Eq.(2.3) then becomes linear in terms of Cj, i = 
1, 2, ..., m. The C's can be determined from the first m equations, or determined 
approximately by applying the least squares method to the encire set of equations. The 
nonlinearity associated in finding e^ a i, which is related to the frequencies joa and -jcOj, is 
concentrated in Eq.(2.4). The above procedure is known as Prony's method. 

Now, in the case of Eq.(2.1), since the roots of Eq.(2.4) occur in reciprocal pairs, 
Eq.(2.4) must be invariant under the substitution of e'J a i for eJ a i, so that we must have 
a 2n = 1( a (2n-l) = a l’ -» a n+ l = a n-l- Thus e T( 2-4) becomes 



ei 2nw + a t e^ 213 ' 1 ) 03 + ... + a n4 e i( n+1 )<° + a n e i n03 
+ a n . 1 ei ^- 1 ) 03 + ... + a { > + 1 =0 
or 

e jnm [ (g jnco + g-jncO) + ai (e j(n-l)co + e -j(n-l)a>) + _ + 

« n .l (el«> + e-i 03 ) + « n ] = 0 

since e^ n03 # 0, we have 



2 cosnco + 2 cos(n-l)o) + ... + 2 a n _j costo + a n = 0 (2.6) 

Now noting that m=2n, and applying the above symmetry of a's to Eq.^2.5) yields 

[Ref. ii] 
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{x(0) + x(2n)} + (x(l) + x(2n-l)} ctj 4- ... + (x(n-l) +• x(n+l)} a n _j 
+ x(n) a n = 0 

{x(l) + x(2n+l)} + (x(2) 4- x(2n)} aj 4- ... + {x(n) + x(n+2)} a n _j 
4-x(n+l)o n = 0 



(2.7) 



(x(N-2n-l) + x(N-l)) + (x(N-2n) 4- x(N-2)} a j 4- ... 4- {x(N-n-2) 
4- x(N-n)} a n _^ 4- x(N-n-l) a n = 0 



Eq.(2.7) consists of a set of N-2n equations and n unknowns. This set can be solved 
exactly for the a's if N = 3n, or solved approximately by the least squares method if 
N> 3n, and then the to's are determined from Eq.(2.6). 

C. ESTIMATION OF FREQUENCY, AMPLITUDE AND PHASE 

If a process under measurement contains an unknown number of sinusoids of 
unknown frequencies and amplitudes, a variant of Prony's method can be used to 
determine the number of sinusoids and their associated frequencies and amplitudes. As 
noted above. Prony's method is a technique for modeling data of equally spaced 
samples by a linear combination of exponentials. An exponential curve having M 
exponentials terms can be determined from the 2M data measurements. Each 
exponential term A^i^ has two parameters — an amplitude A^ and an exponent a ; 
where a. can be real or imaginary. For the case where only an approximate fit w r ith M 
exponentials to a data set ofN samples is desired, such that N>2M, a least squares 
estimation procedure is used. 
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The model assumed is a set of M exponentials of arbitrary amplitude, phase, 
frequency and damping factor. A process consisting of M real undamped (a is 
imaginary) sinusoids can be expressed as 



x(k) - L (az., k + a, V k ) (2-8) 

M 

= £ A. cos(2itfkT + G.) 
i=l ' 1 1 

with a = (Aj/2)e$i and z. = ei 27T ^, where A { is the amplitude, f } is the frequency and 0 ; 

is the phase of the i m sinusoid, respectively, and T is the sampling interval. 

Finding {A;, Gj, fj} and M that minimize the squared error is a difficult nonlinear 

least squares problem. Prony's method solves two sequential sets of linear equations 

with an intermediate polynomial rooting step that concentrates the nonlinearity of the 

problem in the polynomial rooting procedure (similar to Eq.(2.4)). 

* 

Define the polynomial, B(z), which has Z; and z ; as its roots, given by 



M ft 2M -y\A : 

B(z) = II (z - z.)(z - z. ) = £ b. z 2M *I = 0 (2.9) 

i= 1 * 1 

with b Q = 1 and the b- being real coefficients. Since the roots are of unit modulus occur 
in complex conjugate pairs, then Eq.(2.9) must be invariant under the substitution z 
for z. Therefore, Eq.(2.9) can be written as 



z 2M B(l/z) = z 2M 2^ b. zi* 2M = A £ l b. zi = 0 



2M 



[ = 0 i 



1 = 0 * 



( 2 . 10 ) 



Comparing Eqs.(2.9) and (2.10), we may conclude that bj = b 2 ]yj _ j for j = 0,1, ... ,M, 
with bo = b 2 ^ = 1. Thus the requirement for complex conjugate root pairs of unit 
modulus is implemented by constraining the polynomial coefficients to be symmetric 
about the center element. Hence the realization of B(z) is a linear phase FIR filter. 
Based on order 2iM, a linear prediction error can be written as 
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( 2 . 11 ) 



e(k) = L b. [x(k + j) + x(2M-k + j)] 
1 = 0 * 



j=o 



which reduces the number of coefficients required by one-half. 

The coefficients b,. b-,, ... . b_ are determined in a least sauares fashion bv 

i - m 

minimizing the total squared error. 

N-2M o 

E = Z £-(k) (2.12) 

i=0 

which yields the normal equations 
m N-2M 



This equation can be solved recursively for increasing order M. 

From the estimated (b^} values, the {zj are determined using Eq.(2.9). This gives 
the frequency estimates. 




(2.13) 



for k = 1, ... ,M 



f = tan* 1 [Im(z i )/Re(z i )] / 2irT 



(2.14) 



To obtain the (a) a second set of normal equations is solved, 





(2.15) 



for k = 0,l, ... ,M 



17 



The set {a} provides both amplitude (Aj),or power, and phase (0 ; ) information: 



Aj = | a. | (2.16) 

6j = tan' 1 ( Im(b.)/Re(b.) ] (2.17) 

In the foregoing we have shown that the frequency estimation as part of Prony's 
method requires a polynomial with a linear phase property. In the next chapter, we 
show that a lattice structure can be utilized to implement a linear phase transfer 
function. 
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III. LINEAR PHASE FIR FILTERING 



A. INTRODUCTION 

Every finite-impulse response (FIR) filter has two distinct properties [Ref. 5] ; 
first, it is always stable; second, if it is not causal, it can always be made to be causal 
by introducing a finite delay. 

FIR filters can be designed so that their frequency responses have an exactly 
linear phase characteristic. FIR filters with linear phase are important in applications 
like speech processing and data transmission, because in these applications a nonlinear 
phase filter is harmful. 

The impulse response sequence of a linear phase FIR filter exhibits a kind of 
symmetry, e.g., h(n) = h(N-l-n) for n = 0, 1, ..., (N/2)-l (assuming that N is even). 
In general, FIR filters require a large number of coefficients to adequately meet with 
the required filter specifications. The symmetry property of the linear phase FIR 
filters, however, helps reduce the number of coefficients by nearly one half resulting in 
substantial computational savings. 

The various filter realizations, or structures that are frequently considered are the 
direct form, the cascade form, the parallel form, and the lattice form. The lattice form 
realization is of particular interest because of its superior numerical performance, and 
modularity in the structure. Lattice realizations have been successfully utilized in 
filtering and spectral analysis, and in modeling of some physical process like speech, 
geophysical data etc. [Ref. 8], The operation of a lattice filter is completely described 
by specifying the sequence of reflection coefficients that characterize the individual 
stages of the filter. 

In this chapter, we will discuss the basic concepts of the linear phase FIR filter 
and the related lattice structure. And finally, we will show three examples of obtaining 

the lattice structures realizing for both phase and non-linear phase FIR. transfer 
functions (Appendix A). 

B. LINEAR PHASE FIR FILTERS [REF. 3] 

Let (h(n)} be a causal finite duration sequence defined over the interval 
O^n^N-1 having the linear phase symmetry property given by h(n) = h(N-l-n). The 
z transform of (h(n)} can be written as 
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N-l 

H(z) = E h(n)z' n 
n = 0 



(3.1) 



If N is even, then Eq.(3.1) can be rewritten as 



(N/2>1 _ N-l 

H(z) = £ h(n)z‘ n + £ h(n)z* n 

n=0 n=N/2 



Now applying the symmetry property of h(n) in the second term on the right hand side 
yields 



W h(„)z-n 
n=0 

which can be simplified to 



H(z) = {NZ f l h(n)z‘ n hfN-l-n^N- 1 -”) 

n=0 n=0 



H(z) = ^ E^ h(n) {z‘ n + z -(N-l-n)} 
n=0 



(3.2) 



If N is odd, then Eq.(3.1) becomes 



H(z) = t(N J£ H h(n){z' n + z-(N'l-n)] + h(h±) z W-l)/2] 



(3.3) 



If we evaluate Eqs.(3.2) and (3.3) for z=el (0 , we obtain the discrete Fourier transform 
of the filter sequence h(n) defined as 



= N £ A h(n) e'i 0511 
n = 0 



(3.4) 



For the case when N is even, we then have the discrete Fourier transform 
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(3.5) 






e -jco{(N-l)/2} 



(N/2)-l 

L 2h(n)cos( to(n-(N-l)/2)}J) 
n = 0 



and for N odd, we obtain 



liCeJ 03 ) ~ [h {( W- i )/ 2 } + (N f 3)/2 2h(n) cos[oj(n-(N-l)/2) ]]. 

n = 0 



(3.6) 



In both cases above, the sums in brackets are real, implying a linear phase shift 
corresponding to a delay of ( N- 1 )/ 2 samples. Figures 3.1 and 3.2 show the direct form 
realizations of an FIR filter with linear phase for both N even and N odd cases. 




I 



Figure 3.1 Direct-form realization for an FIR filter with linear phase (N = even). 
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Figure 3.2 Direct-form realization for an FIR filter with linear phase (N= odd). 

C. LATTICE FILTERS [REF. 4[ 

The basic single section lattice structure is shown in Figure 3.3 where x(n) is the 
x(n) and fj(n) and gj(n) are generally known as the forward and backward prediction 
errors, respectively. The defining equations of the lattice are given by 

f 0 ( n ) = £o( n ) = x ( n ) 

f,(n) = f 0 (n) + k lgo (n-l) (3.7) 

g t (n) = k^n) - g 0 (n- 1) 

where kj is called the lattice reflection coefficient. If another section is cascaded with 
the first one the resulting equations for the next order prediction errors are given by 
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(3.8) 



f 2 (n) = fj(n) + k 2 gj(n-l) 
g 2 (n) = k 2 f 1 (n) + gjCn-l) 

Substituting for fj(n) and gj(n-l) from Eq.(3.7) yields 

f>(n) = x(n) + (kj + ^k-Jxfn-l) +■ k^x(n-2) 
or 

f 2 (n) = x(n) + b 2 ix(n-l) + b 2 2 *(n- 2 ) 



and 

g 2 (n)= k 2 x(n)+ (kj + kjk^xCn-lJ + xOi^) (3.9) 

or 

g 2 (n) = b 2)2 x(n) + b 2)1 x(n-l) + x(n-2) 
where b 2 i = kj(l + k 2 ) and b 2 2 = ^2* 

Thus, by induction, for a cascade connection of M lattice sections we have the general 
expressions 



M 

f M< n ) “jJ-q b M,i x ^ n * 1 ) 

M 

=4 b M,M-i x ( n_i ) 



(3.10) 



tVi 

where b^,j q= 1. Eq.(3.10) represents the FIR filter type equations to obtain the M n 
order forward and backward prediction errors. Now taking the z-transform of 
Eq.(3.10) yields 
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(3.11) 



F M< z > = ^11) b M,i z '* X(z > 
C M (z) = .]~ Q b M ( M-i z l X(z) 



(3.12) 



For an unit impulse input, i.e., when X(z)= 1, F^(z) and G^j(z) represent the forward 
and backward prediction error filter transfer functions, respectively, and 



G m (z> = z' m F m (z-') 



(3.13) 




Figure 3.3 Lattice Section. 

By inspection of Figure 3.3 we see that the M tb stage lattice reflection coefficient 
is equal to the FIR lilter coefficient b^j ^ in Eqs.(3.11) and (3.12). In other words, 
we have 
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k M“ b M,M 



(3.14) 



In fact, the FIR like prediction error filter coefficients can be iteratively calculated 
starting from Eq.(3.14). The algorithm, to calculate the filter coefficients •, also 
makes use of the well known lattice property that an order lattice contains all 
prediction filters of orders m^M. Now consider the m tb section (l<m^M) in the 
cascade connection of M lattice sections which can be described by the following 



equations: 

F m^ z ) = F m-l( z ) + k m z ’ lG m-l( z ) (3-15) 

G m( z ) = k m F m-l( z ) + z ' lG m-l( z ) ( 3 - 16 ) 

Eq.(3.I6) can be rewritten as 

G m-l( z ) = zG m ( z ) ' zk m F m-l^ z ) (3-17) 

Substituting Eq.(3.17) into Eq.(3.15) yields 

F m^ z ^ — F m-l^' l ” k m^ G m^’ k m F m-l^ z ^ (3.18) 



Therefore, the (m-l) m order forward prediction error transfer function can be written 

tin 

in terms of the nr 11 order forward and backward prediction error transfer functions as 
follows: 



F m-l( z ) = 



F m( z > ~ k m G m( z > 



1 - k' 



(3.19) 



m 



where k m =b mrn ?“ 1. Recalling Eq.(3.13), that is, 



G m( z ) = z‘ m F^z' 1 ) 



rm 



(3.20) 



By substituting Eq.(3.20) into Eq.(3.l9) we find that 
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F m-l( z ) 



(3.21) 



F m ( z ) * k 



m 



z m F m ( z' 1 ) 



1 - k" 



m 



Thus, from Eq.(3.21) we see that the next lower order polynomial F m _j(z) can be 
calculated knowing F m (z). Following the foregoing procedure we can find 
k m _j = b m _j from F m _j(z), and also obtain F m _ 2 (z)- This way, for a given M tF 
order polynomial Fjyj(z), we can determine all the lattice reflection coefficients k m , 
m=M, M-l, ..., 2, 1. This procedure is known as the step-down procedure [Ref. 1], 
and can be summarized as follows: Let the given M tF order polynomial be 



f m( z > = 




b 



M,i z 



-l 



(3.22) 



and by replacing M with m we have an m 



th 



order the polynomial for m lattice sections 




(3.23) 



where m= M, M-l, ..., 2, 1. We now define another polynomial given by 



F n/ Z ) ^m,i z 



i 



(3.24) 



As we step through the procedure from m= M to m= 1 Eq.(3.24) can be expressed as 



- (z' 1 ) = f 1 b -z 111 * 1 
nr 1 -=-q u m.m-r 



(3.25) 



Define k m = b mm and obtain the m-l lF order polynomial as follows: 
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F m-l( z ) = 



(3.26) 



F m( z ) ' ^m z F m ( z ) 



1 - k" 



m 



Substituting Eqs.(3.23) and (3.24) into Eq.(3.26) yields 



w 1 - 



,J5 



b mi z "‘ * k m z '“'^ ) b m,m-i zra '‘ 



m 



^ 0 b m-l,i z 1 



^b— :Z*i - k,^ m ;Z"* 

• x m,i m j m,m-i 



m 



(3.27) 



Then, by equating the coefficients of like powers of z on both sides of Eq.(3.27), we 

• + U 

obtain the computational expression for the (m-1 r n order polynomial coefficients as 



b m-l,i 



b mi ' ^m b 



m.m-i 




(3.28) 



where m = M, M-1, ..., 1 and i = 0, 1, ..., m-1, 1^ = b mm and |k m | < 1 for a 
minimum phase polynomial F m (z). 

D. LINEAR PHASE FIR LATTICE FILTER 

In the foregoing we considered obtaining a lattice structure from a given FIR 
transfer function, and vice versa. In what follows we will deal with a special case of 
FIR filtering, namely, the linear phase FIR filter. Let (h(n)} be a causal finite duration 
linear phase sequence defined over the interval 0<n^N-l. From Eqs.(3.2) and (3.3), 

the z transform of (h(n)} can be written as 



H N-l( z ) = f M< z ) + Z ‘ D g m( z ) 



(3.29) 
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Figure 3.4 Linear Phase FIR Lattice Filter. 

where F-^(z) and Gj^j(z) are forward and backward filter transfer functions, 
respectively, and D represents the number of unit delays. We now consider the N odd 
and N even cases separately. 

iV odd case : For N odd, Eq.(3.29) can be written as 



H N-l( z > = a 0 + a l z * 1 + a 2 z * 2 + ••• + a (N-3)/2 z 



-(N*3)/2 , _ _-(N-l)/2 

+ a (N-l)/2 z 



+ a (N-3)/2 z 



■(N-D/: 



a 2 z-( N - 3 > + a ( z-( N - 2 ) - 



= a 0 + ajz* 1 + a 2 z' 2 + ... + a ( N-3)/2 z * (N * 3)/2 + 2 (1/2) a (N . 1)/2 z' (N ' 1)/2 
4- a (N . 3)/2 z-(N+l)/2 + ... + a2Z -(N-3) + ^-(N.2) + a()Z -(N-l) 
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Note that we are splitting the coefficient 
(l/ 2 > a (N-l)/2 each - 



into two coefficients of value 



H N . j ( z ) = a Q + a t z _1 + a 2 z' 2 + ... + a ^ N . 3 ) /2 z ^ N ‘ 3 ^ 2 (3-30) 

+ (1/2) a ( N . 1)/2 z -( N * 1 )/ 2 + z *( N * 1 )/2 [(i/ 2 )a (N _ 1)/2 

+ a ( N . 3 )/ 2 z ' 1 + ••• + a 2 z *^* a ^ 2 4 - ajz *^* 2 )/ 2 
4- a 0 z -( N - 1 >/2 ] 

H n .!(z) = F m (z) + z-( n - 1 )/2 Gm(Z ) (3.31) 

where the order of the polynomials F^(z) and G^fz), M = (N-l)/2, and the 
number of unit delays, D = (N-l)/2, the forward transfer function, 



F m ( z ) = a 0 + a lZ * 1 4 - ... + a ( N . 3 )/ 2 z -( N - 3 )/ 2 + ( l ^ a ^. j ^ z ^- 1 )/ 2 



Fm^ 2 ' 1 ) “ a 0 + a l z + - + a (N-3)/2 z ^ N " 3 ^ 2 + ( 1 / 2 ) a (N-l)/2 z( ' N ’ 1 ^ 2 



and the backward transfer function, 



G m ( z ) = z-(N-0/2 F m (z‘ 1 ) = (l/2)a (N . 1)/2 - a (N . 3 )/2 z' 1 



(3.3-!) 



+ 



4 - ai z -( N - 3 >/ 2 - Haoz - CN * 1 )/ 2 



(3.32) 



(3.33) 
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N even case : For N even, Eq.(3.29) can be written as 






H n .j(z) 



= a 0 + ajz' 1 + a 2 z* 2 + ... + a (N _ 2)/2 z-( N - 2 >/ 2 

+ ... + a 2 z'<*-3) + a lZ -( N - 2 > - aoz'^* 1 ) 

= a Q + ajz* 1 + a 2 z* 2 + ... 4- a(N- 2 )/ 2 z ' (N ' 2)/2 

^• N/2 [ a (N.2V2^-^^ 

+ a 1 z -(^)/ 2 + a 0 z ^- 2 )/ 2 ] 



+ a (N-2)/2 z 



-N/2 



(3.35) 



H N-1 (z) = F m (z) + z' N ' /2 G m (z) ( 3-36) 

where the order of the polynomials F-yj(z) and G^(z), M = {(N/2)-l}, and the 
number of unit delays, D = N/2, the forward transfer function, 



f M ( z ) = a 0 + a l z4 + a 2 z " 2 + - + a (N-2)/2 z ' (N * 2)/2 ( 3 - 37 ) 

f m ( z_1 ) = a 0 + a l z + a 2 z2 + - + a ( N -2)/2 z ^ N * 2 ^ 2 ( 3 - 38 ) 

and the backward transfer function. 

g m (z) = z-( N - 2 V 2 F m (z' 1 ) = a (N _ 2 y 2 + ... + a 2 z*( N * 6 ) /2 (3.39) 

4 - ... 4 - a lZ -( N - 4 )/ 2 + a 0 z -( N - 2 )/ 2 
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Figure 3.4 shows a linear phase FIR lattice realization where the filter output is 
obtained by adding the order forward prediction error, and the order delayed 
(by D unit delays) backward prediction error. Combining the discussion in this section 
and that in Section (III.C), we can summarize the algorithm for converting a given 
FIR transfer function into a lattice structure as follows: 

(i) Find if N is even, or odd 

(ii) For N odd: M = (N-l)/2, D = (N-l)/2 

(iii) For N even: M = {(N/2)-l}, D = N/2 

(iv) From Fjyj(z), obtain M reflection coefficients as discussed in Eqs.(3.22) to (3.26) 

(v) Implement the lattice as shown in Figure 3.4 

E. LATTICE REALIZATION OF A GENERAL FIR TRANSFER FUNCTION 

In the foregoing we have considered a polynomial Fjyj(z) with b^ q= 1. However, in 
general we have b^o^L in this section, we shall modify the lattice realization 
algorithm presented in Section (III.C) to suit an arbitrary FIR transfer function with 
0 ^ * [Ref- 22]. The m^ order polynomials F m (z) and G m (z) can be obtained 
from Eqs.(3.15) and (3.16): 

F m< z > _ s m F m-l< z > + k m z '' G m-l( z ) < 3 '40) 

G m< z > “ k m F m-l( z > + s m z '' G m-l( z ) < 3 - 41 > 

where the coefficients s m = b m q and k m = b m m are recognized as the reflection 
coefficients. Eq.(3.41) can be rewritten as 

G m-l< z > “ s '‘m z G m< z > ' s 'm k m z F m-l( z > < 3 ' 42 > 

Substituting Eq.(3.42) into Eq.(3.40) yields 

F m< z > ' s m F m-l( z > + s '‘m k mf G ,J z > ' koiPm-lf 2 '* 

Therefore, the (m-i) 1 ^ 1 order forward prediction error 
written in terms of the m^ order forward and backward 
functions as follows: 



(3.43) 

transfer function can be 
prediction error transfer 
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Figure 3.5 FIR Lattice. 






— “m m 



F_(z) - k m G m( z ) 



s 2 - k l 
m m 



(3.44) 



where s m ^k m . Substituting Eqs.(3.20), (3.23) and (3.25) into Eq.(3.44) yields 




-i _ k 7 -m^ b m-i 
'jn^j D m,i z K m z j^pn,m-i z 



1-2 



m 



m 




m i i 0 b m,i z 



' b m,m-i z 




(3.45) 
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From Eq.(3.45), we obtain the computational expression for the (m-l) 1 ^ order 
polynomial coefficients is 



5 m- 1 ,i 



s m^m.i ' ^m^m.m-i 



m 



m 



( 3 . 46 ) 



where m = M, M-l, .... 1 and i = 0, 1, .... m-1, k m = b mm , s m = b m0 and s m > k m 
for a minimum phase polynomial F m (z). It may be noted that s m = i for m = 1, 2, 
M-l. This indicates that we have only one s-coefficient, i.e., s^, which requires an 
extra multiplication in the M in lattice section as shown in Figure 3.5. 
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IV. LMS ALGORITHM 



A. INTRODUCTION . 

Conventionally an adaptive filter is composed of a tapped delay line or 
transversal structure with adjustable coefficients or weights and an adaptive algorithm 
which updates the coefficients continuously based on some performance criterion. 

The design of a fixed coefficient filter is based on the prior knowledge of both 
signal and noise. Adaptive filters, on the other hand, have the ability to adjust their 
own parameters automatically, and their design requires little, or no a priori knowledge 
of signal or noise characteristics (Ref. 14]. However, the designer has to choose the 
order of the filter and the type of the algorithm. Also the adaptive filter usually 
requires a large initial transient time (i.e., the initial filter convergence period). 

For stationary stochastic inputs, the mean square error, the difference between 
the filter output and an externally supplied input called the "desired response", is a 
quadratic function of the filter coefficients, a paraboloid with a single fixed minimum 
point that can be sought by gradient techniques (Ref. 16]. 

In the previous chapter we showed that the operation of a multistage lattice filter 
is completely described by specifying the sequence of reflection coefficients that 
characterize the individual stages of the filter. In this chapter we briefly discuss the 
least-mean-square (LMS) adaptive algorithm that results from attempting to minimize 
the mean square error and present a summary of some LMS algorithms for lattice that 
have been reported in the literature. Also included are the computer simulation results. 

The least-mean-square (LMS) adaptive algorithm minimizes the mean square 
error £(k) by recursively altering the filter coefficient vector B(k) at each sampling 
instant. The original Widrow-Hoff LMS algorithm is B(k+ l) = B(k) + 2fie(k)X(k), 
where n is a convergence factor controlling stability and rate of adaptation (Ref. 15]. 
The algorithm is based on the method of steepest descent, moving b(k) in proportion 
to the instantaneous gradient estimate of the mean square error. 

When the filter input is stationary, the backward prediction errors are orthogonal 
to each other, with the result that the successive stages of the lattice filter are 
decoupled from each other [Ref 8]. This means that the global optimization of a 
multistage lattice filter may indeed be accomplished as a sequence of local optimization 
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problems, one at each stage of the lattice filter. Accordingly, it is a straightforward 
matter to increase the order of the lattice filter by simply adding one or more stages 
without affecting the earlier design computations. The successive orthogonalization 
provided by the lattice offers advantages in adaptive convergence rate which cannot be 
achieved with tapped delay lines [Ref. 17,18]. 

B. SUMMARY OF THE LMS ALGORITHM 

The LMS algorithm uses an estimate of the gradient of the mean square error obtained 
from the adaptive linear combiner which is a combination of a transversal structure 
and an adder. The adaptive linear combiner can be shown in two basic ways, 
depending on whether the input is available in parallel form (multiple inputs), or in 
serial form (single input) [Ref. 6]. 
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In the following we present a brief derivation of the LMS algorithm. Let us choose the 
single input form, then the filter output is given by 

y(k) = b 0 (k)x(k) + b 1 (k)x(k- 1) + ... + b N . 1 (k)x(k-N+ 1) 

N-l 

= £ bj(k) x(k-i) 

i = 0 

= X T (k) B(k) 

= 5 T (k)X(k) (4.1) 

where X and B are the input signal vector and the filter coefficient vector, respectively 
and (N-l) is the order of the filter. The input signal vector X and the filter coefficient 
vector B arc defined as 
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m = 





m = 


m ' 


x(k-N+ 1) 




^N-l« 



Therefore, the output y(k) is equal to the inner product of X(k) and S( k). 

The error e(k) is defined as the difference between the desired response d(k) and 
the actual response v(k), 



e(k) = d(k) - X T (k)S(k) = d(k) - S T (k)X(k) (4.2) 

The purpose of the adaptive algorithm is to adjust the filter coefficients of the adaptive 
linear combiner to minimize the mean-square error. A general expression for mean 
square error as a function of the filter coefficient values, assuming that the input 
signals and the desired response are statistically stationary and that the filter 
coefficients are fixed, can be derived in the following manner [Ref 6]. The squared 
error is 



e 2 (k) = d 2 (k) - 2d(k)X T (k)S(k) + B T (k)X(k)X T (k)S(k) (4.3) 

Taking the expected value of both sides yields the mean square error, 



E[e 2 (k)] = E[d 2 (k)] - 2 E[d(k)X T (k)] B(k) + B T (k) E[X(k)X T (k)] B (4.4) 

Defining the vector P as the cross-correlation between the desired response and the 
input vector, we have 



P = E [ d(k) X(k) ] (4.5) 

Similarly the input autocorrelation matrix R is defined as 
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(4.6) 



R = E [ X(k) X T (k) ] 

Thus the mean-square error can be expressed as 

£(k) = E[e 2 (k)| = E[d 2 (k)] - 2 P T 5(k) .+ B T (k) R B{ k) (4.7) 

The gradient Vs(k) of the mean square error function is obtained by differentiating 
Eq.(4.7) with respect to the filter coefficient vector as follows: 



V£(k) = 



dE[e 2 (k)] 

5b 0 (k) 

<3E[e 2 (k)] 

3b N .,(k) 



= -2 P + 2 R 5(k) 



(4.3) 



The LMS algorithm is an implementation of the method of steepest descent. 
According to this method, the next filter coefficient vector is equal to the present filter 
coefficient vector B( k) plus a change proportional to negative of the gradient, V£(k): 



S(k+ 1) = 5(k) - V£(k) (4.9) 

The parameter p is the factor that controls stability and rate of convergence. In other 
words, the first term on the right hand side consists of the past information and the 
second term represents the new, or updated information. 

The LMS algorithm estimates an instantaneous gradient in a crude but efficient 

manner by assuming that e~(k), the square of a single error sample, is an estimate of 
the mean-square error and by differentiating e 2 (k) with respect to 3(k). The estimated 

gradient is given by the following expression 
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(4.10) 



7c(k) = 



<5e 2 (k) 




<5e(k) 


ob 0 (k) 




db 0 (k) 




= 2 e(k) 


• 


3e 2 (k) 




de(k) 






The estimated gradient components are related to the partial derivatives of the 
instantaneous error with respect to the filter coefficient components. Thus the 
expression for the gradient estimate can be simplified to 



Vs(k) = -2 e(k) X(k) (4.11) 

Using this estimate in place of the true gradient in Eq.(4.11) yields the Widrow-Hoff 
LMS algorithm 



g(k+l) = S(k) + 2 n e(k) X(k) (4.12) 

Since the filter coefficient changes at each iteration are based on imperfect gradient 
estimates, we would expect the adaptive process to be noisy, that is, it would not 
follow the true line of steepest descent on the performance surface. The LMS algorithm 
can be implemented in a practical system without squaring, averaging, or 
differentiation and is elegant in its simplicity and efficiency. Each component of the 
gradient vector is obtained from a single data sample without perturbing the filter 
coefficient vector. 

C. ADAPTIVE LATTICE ALGORITHMS 

The parameters to update in a multistage lattice are its refection coefficients. 
Several algorithms have been proposed in the literature for updating the lattice 

reflection coefficients [Ref. 7,8,17,18,20,23-28]. In this section we briefly summarize 
some of those algorithms. Consider a lattice filter of order M. For stage m of the 
filter, the flow of signals is described by the following pair of equations. 
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Figure 4.3 A Single Stage of Lattice Filter. 

« k > - + (4 , 3) 

Sm< k >- em-l< k -'> + k m,b f m.l< k ) 

where m= 1,2, ... ,M 

For stage m, the forward reflection coefficient is denoted by k m ^ and the backward 
reflection coefficient is denoted by k m f m (k) and g m (k) are the forward prediction 
error and backward prediction error of stage m respectively. Before presenting the 

adaptive algorithms for the lattice, let us quickly summarize some nori-adaptive 
reflection coefficient estimation methods. Later on we can obtain the adaptive 
updating equations as approximations to these methods. 
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1. Non-adaptive Methods 

When the lattice coefficients have fixed, non-adaptive values, several methods 
have been proposed for computing these values as functions of the correlation statistics 
of the input. Two of these which are based on mean-square error (MSE) minimization 
[Ref.28] are given below. 

Method 1 : In this method, we choose k m ^ to minimize the mean forward 

prediction error power, E[f 2 m (k)] and ^ to minimize the mean backward prediction 
error power, E[g 2 m (k)]. Here we give the final equations for k m ^ and k m ^ without 
going into the details. The forward and backward reflection coefficients are given by 



_ Elf^Wg^tk-l)] 
m ’ f E[g 2 m-1 (k-l)] 

Elf m .,(k)g m .,(k.l)l 
m ' b E[f 2 m .,(k)l 



(4.14) 



where both k m p and k m ^ are obtained from the cross-correlations between the 

fU 9 9 

(m-l) in order forward and backward prediction errors normalized by respective 
prediction error powers. 



Method 2 : For a single channel lattice using real data and coefficients we can, 
however, show that k m p=k m ^=k m . Based on the condition that k m f=k m ^ = k m , 
we can now minimize either E[f 2 m (k)], or E[g 2 m (k)] in order to obtain an optimum k m . 
However, it seems more logical to minimize the sum E[f 2 m (k)] + E[g 2 m (k)] as suggested 
in [Ref. 28]. The resulting reflection coefficient equation in its final form is given by 



2 Etf^wg^ck-ni 

^^m-l^)]"^ E[g 2 m . 1 (k-l)] 



(4.1b) 



where we have normalized the cross-correlation term with respect to the sum of mean 
prediction error powers. In both of these methods, the coefficients at stage m can be 
computed independently of those following that stage. Thus, optimum values can be 
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computed successively along the structure without affecting the other coefficients. This 
phenomenon lends to the modular nature of the lattice structure. 

2. Adaptive Methods 

An instantaneous, gradient-descent adaptive algorithm minimizes a mean- 
square error criterion can be derived for a generalized problem. An adaptive lattice 
structure is suggested by incorporating time varying coefficients k m ^k) and k m ^(k) in 
the lattice and generating algorithms for the two methods described above. The 
resulting procedures are : 

Method I : Corresponding to method 1 of non-adaptive procedures metioned 
above, we can obtain the following update equations 



W k + = W k > + 






m- 



i,f00 



t f m( k > Sm-l^ 1 )] 



K m,b( k ' ^ ^m-1® 

c m-l,gW 



(4.16) 



where |i is an adaptive step size parameter, X is a positive weighting constant, and the 
forward power estimate at the (m-l) tk stage, f(k), is 

* 2 m.l,(< k ) ” k °W k -D + d-J-) ts 2 m-lC k -DJ 

and the backward power estimate at the (m-l) tk stage, > 0 (k), is 



Method 2 : Eq.(4.15) can be approximated to obtain a recursive update equation 

as follows: 



lhn(k+ 1) = k m (k) + m [f m (k)g m .,(k.l) + r m .,(k)g m (k)] (4.17) 
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For the basic structure of the single-channel lattice, ie, k m p = k m ^ = k m , the reflection 
coefficients or filter coefficients k m (k) may be updated using a modification of the 
LMS algorithm of the LMS algorithm [Ref. 17,20]. 



where X is a positive weighting constant satisfying the criterion 0<X^ 1 then controls 
the bandwidth of the filter and the resulting power averaging time. A power estimate 
is required at each stage in the lattice due to the fact that the forward and reverse error 
sequences have decreased power with increased stage number. 

Method 3 : A third successful method of implementing an adaptive algorithm for 
FIR lattice structures has been reported by Griffiths [Ref. 17,18]. This algorithm has 
been originally discussed for a noise canceller application. The form of the lattice 
noise-cancelling adaptive filter is shown in Figure 4.4. The lattice noise-cancelling 
adaptive filter consists of an M stage linear prediction lattice for the reference signal 
x r (k) together with a set of tap coefficients V m (k) which provide the noise-cancelling 
subtraction paths. Griffiths' algorithm is briefly presented in the following. The 
update equations for k m (k) are given by 




where ^“m.^k) is the power estimate at the (m-l)^ stage. Now the updated power 



estimate is 



<7 2 m .i(k) = X <r 2 m .i(k-l) + (1-X) [g^fk-l) + f^k)] 



(4.19) 




where the power estimate at the (m-l) 1 * 1 stage, <7 2 m _j(k) is 



(T^fk) = Xtr^k-l) + (1-X) [g^k-l) + f 2 m .!(k)] 



(4.21) 
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Figure 4.4 Lattice Form Implementation of Noise-cancelling Filter. 

Each coefiicient V n (k) can be determined independently of V n (k) for n> m, because of 
orthogonalization provided by the lattice. Thus the resulting algorithm is 



v m< k + » - v ,„( k > + T /- k ' l‘m< k > 8m( k )l (« 2 > 

T m \ K ) 

where associated power measurement y 2 m (k) is 

Y 2 m (k) = Xy 2 m (k-1) + (1-X) (g m (k-l) g m (k-l)] (4.23) 

and £ m (k) is the stage error signal as shown in Figure 4.4. 
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D. SIMULATION RESULTS 




Figure 4.5 System Identification Experiment. 

For the purpose of computer simulation, let us consider an approximate 
algorithm given by 



k m< k+i > “ k m< k ) + „ z 1 *,,' lg m -i( k -l)l e(k) 



(4.24) 



where 



<r 2 m (k) = X<r 2 m (k-1) + (1-X) gv^k-l) 



(4.25) 
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The configuration used in the simulation is a system identification experiment as shown 
in Figure 4.5. The fixed filter transfer function considered is given by 
H(z) = 1 - 0.89 z' 1 + 0.25 z 2 

The convergence performance of the LMS algorithm (as given by Eq.(4.24)) can be 
observed by plotting the error, e(k). versus the update iteration, k. called learning 
curves. The input x(k) to both fixed and adaptive filters is a white noise sequence with 
zero mean and unit variance. Figures 4.6 to 4.8 show the learning curves for the above 
example where we have used three different values for the adaptation constant, ji. In 
the next chapter, we derive a new algorithm for updating the reflection coefficients, 
based on the least mean square principle and the steepest descent algorithm. 
Improvement in convergence speed will be shown using the new algorithm. 




Figure 4.6 Learning Curve (ft = 0.01). 
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Figure 4.7 Learning Curve (jj. = 0.1). 




Figure 4.8 Learning Curve (fi = 0.5). 
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V. ADAPTIVE LINEAR PHASE LATTICE ALGORITHM 

A. INTRODUCTION 

In Chapter III we dealt with the realization of fixed coefficient FIR lattice filter with 
both linear and nonlinear phase characteristics. We reviewed the basics of LMS 
algorithm and some adaptive lattice realizations reported in the literature in the 
previous chapter. The three adaptive lattice algorithms discussed are direct 
approximations of their non-adaptive counterparts. None of them estimates a gradient 
as required by the LMS algorithm. In this chapter we present a derivation for a new 
LMS based FIR adaptive lattice algorithm and extend it to the linear phase case as 
well. The new algorithm is then used in the estimation of spectral lines in white noise. 
Results of simulation are included. 

B. LMS ALGORITHM FOR THE FIR LATTICE 

From Eqs.(3.40) and (3.41) the FIR lattice equations can be written as 



s m + ^m^m-l^'^) 

+ < 5 - 2 > 

where m= 1,2,...,M, s m = 1 for m^M and k m are the lattice reflection coefficients. A 
realization of Eqs.(5.1) and (5.2) is shown in Figure 5.1. The lattice input is 
x(k) = gQ(k) = f()(k). The output of the filter is 



y(k) = f M (k) 

= S M ^M-l^) + gM-l(k*l) 

Therefore, the error, e(k), is given by 



e(k) = d(k) - y(k) 

= d(k) - s M f M .](k) - k M g M .i(k-l) 



(5.4) 
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Figure 5.1 Adaptive FIR Lattice Filter. 

where d(k) is the desired signal. The objective is to minimize the mean square error 



J = E { e 2 (k) } 



(5.5) 



The gradient of J with respect to k m and s m , respectively, are given by 




3v(k) 

-2 E(e(k) } 

Uk; 



(5.6) 



and 
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5J 



= -2 E{e(k) f M .j(k)} 



(5.7) 



5s 



M 



Then the LMS algorithm can be formulated as follows: 



dJ 



kj(k+l) = kj(k)-n k ( ^ 



5y(k) 

= kj(k) + 2 n k E(e(k) ( )} 



(5.8) 



_ 5y(k) 

- k j( k) + 2 Mk e(k) ( ) 



and 



5J 

s M (k+!) = s M (k) - \i s ( — — ) 

ds M 

= s M (k) + 2 E{e(k) f M .j(k)} (5.9) 

= s M (k) + 2 n s e(k) f M .j(k) 

where n k and are the adaptation constant, and we have replaced the true gradient 
by its instantaneous estimate. Defining 



z(k) 



5y(k) 




(5.10) 



yields 



kj(k+ 1) = kj(k) + 2 e(k) z(k) 
s M (k+!) = s M (k) + 2 n s e(k) f M .j(k) 
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where j=l,2,...,M. The next step is to estimate the gradient vector z(k). For this 
purpose, define 







and 



(5.12) 







(5.13) 



Then from Eq.(5.3) z(k) is given by 



z(k) = 0 Mj -(k) 

= s M (k) 0 > M .i j(k) + k M (k) , P M . lf j(k-l) + gM-i(k-l) s M.j 



(5.14) 



where 5^ j = (dk^/dkp. Substituting Eqs.(5.1) and (5.2) to Eqs.(5.12) and (5.13) then 
we have 









dk: 



'J’i.jW 



S i (kH S— + k i (k) iiJ— + f i- 



therefore. 



®y(k) - Sj(k)Oj_j_j(k) + kjWFj.^fk-l) + g^fk-DSy 
'J'ij(k) - s i (k)'F i . lij (k-l) + ! j(k) + f M (k)Sy 



(5.15) 

(5.16) 



where 
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(5.17) 






is the Kronecker delta function, and i= 1,2, ... ,M and j= 1,2, ... ,M. If i=0, then 

Eqs.(5.15) and (5.16) are 



Oo.jW ’ 



dx(k) 

<3k ; 



= 0 



(5.18) 



T oj (k) = <D 0( j(k) = 0 (5.19) 

where j = 1,2, ... ,M. Figure 5.2 shows the computation of gradient elements for the 
adaptive FIR lattice algorithm. 

From the Eqs.(5.l5), (5.16) and (5.19) we have 



®ij00 = Tj j(k) = 0 (5.20) 

where l<i<j-l. The computations of gradient elements at each case of j = M,M-l, ... 
,1 are as follows: 



Case j= M : 



^0,M^ “ “ 0 

= ^i(k)0>i.i M (k) + k i (k)'P i . LM (k-l) + 

+ f i-l( k ) 5 i,M 

where i= 1,2, ... ,M. From Eq.(5.20), every gradient element equals to zero except final 
stage elements and gradient elements of the final stage are 
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Figure 5.2 Computation of Gradient Elements for the Adaptive FIR Lattice Algorithm. 



Applying the Eqs.(5.17) and (5.20) yields 



53 







(5.21) 

(5.22) 



Case j= M-l : 

®0,M-l< k > - '•'o.M-lW " 0 

®i,M-l* " s i( k ) 0 i.l,M-l( k ) + k i( k )' J ’i-l,M-l( k -» + «i.iC k -0»i,M-l 
1'i.M.im - ^^^i-l.M-lCk-O + k i( k )®i.l,M.l( k ) + f i-l< k > 6 ,M-l 

where i= 1,2, ... ,M. Applying the Eq.(5.20), last two stage elements are considerable, 
if i= M-l, then we have 

^M-LM-lW = s M-l (k ^.M-2,M-l® + ^M-I^^M^M-l^’ 1 ) + 

T M-l,M-l( k ) = s M-lW' I 'M-2,M-l( k ’ 1 ) + k M-l( k ) <t> M-2,M-lW + f M-2 (k ^M-l,M-l 
Applying the Eqs.(5. 17) and (5.20) yields 



^M-l.M-l^) “ §M-2( k ' 1 ) ( 5 - 23 ) 

T M-l,M-l (k ) = f M-2 (k ) ( 5<24 ) 

And the last stage terms are 

' D M,M-i( k > = s M( k ^M-l,M-l( k ) ~ k iViW' p M-l,M-f k ' 1 ) ~ 

T M,M-l( k ) = s M^ k ^'^M-l,M-l( k ' 1 ) + k M( k ^M-l,M-l( k ^ + f M-l^M,M-l 
Using the Eqs.(5.17), (5.23) and (5.24) yields 
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~ s M^)§M-2( k ' 1 ) + k M (k ) f M-2( k * 1 ) 
= s MW f M-2^') + k M^ k ) SM-2^* 1 ) 



(5.25) 

(5.26) 



Finally, 



Case j= 1 : 



*0,1^ " ^0,l< k > - 0 




(5.27) 

(5.2S) 



where i= 1,2, ... ,M. 

The LMS algorithm derived in the foregoing can be summarized in Table 1. 

Simulation Results : 

The performance of the algorithm summarized in Table 1 has been observed by 
computer simulation. The configuration used is the system identification experiment as 
discussed in Section (IV. D) using the same fixed coefficient filter. The learning curves 
obtained using the new algorithm are shown in Figures 5.3 to 5.5. Comparing these 
with the learning curves in Figures 4.6 to 4.8, which are obtained using approximate 
algorithms, we observe significantly faster convergence rate for the new algorithm. 
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TABLE 1 

LMS ALGORITHM FOR THE FIR LATTICE 

Initialization: 

K( 0) = 0 

s M ( 0 ) = 1 

Lattice : 

x(k) = f 0 (k) = g 0 (k) 

f m< k > = s m < k > f m-l< k > + V k )9 m -l< k - 1 > 

9m< k ) " s n,< k )9 m -l< k -D + V k ) f m-l< k ) 

m = 1, 2, . . . , M 

y(k) = f M (k) 

Update Equations : 

kj(k+l) = kj(k) + 2 [H k /<* 2 k (k)l e(k) Zj (k) 
s M (k+l) = s M (k) + 2 [ H s /cr 2 s ( k)] e( k) f^^k) 
where ji k and fi g are the adaptation constants, 

<* 2 kj (k) = X <J 2 kj (k-l) + (1-X) <J> 2 M/j (k) 
and 

<y 2 s (k) = X o 2 s ( k-1 ) + (1-X) f 2 M . 1 (k) 

are estimations of power in Zj(k) and fjyj_i(k), 
respectively and X is a positive weighting constant, 
0<X< 1. 

Gradient Vector Elements : 



®0,j< k > 


= * 0/J <10 = 0 




®i,J< k ) 


= s i (k)0> i _ 1/ j(k) + 


k i( k )Ti_i, j( k -D 




+ g i-l( k ~ 1 ) < 5 i / j 






= 3 i (k) , F i . 1/j (k-i) 


+ k i (k)O i . 1/j (k) 




+ f i-i< k ) s i,j 




Zj(k) = 


4 > M/j ( k ) 
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Figure 5.3 Learning Curve (ji = 0.1). 




Figure 5.4 Learning Curve ([i = 0.3). 
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Figure 5.5 Learning Curve (p = 0.5). 



C. LINEAR PHASE FIR LATTICE ALGORITHM 

We now extend the LMS algorithm derived in the previous section to the linear 
phase FIR lattice Filter. 

From Eq.(3.29), output of the linear phase FIR lattice can be written as 



y(k) = f M (k) + g M (k-D) (5.29) 

The lattice input is x(k)= fQ(k)= gQ(k). Substituting Eqs.(5.1) and (5.2) in the output 
equation of the Filter yields 



y(k) - s M f^lk) - k M g M _!(k-l) + s M g M -i(k-D-l) + k M f M-l^ k ’ D L 30) 
= S M t f M-l^ k ) + §M-l(k-D-l)] + k M i f M-l( k ' D ) + 

and the error, e(k), is given by 
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d(k) 




Figure 5.6 Adaptive Linear Phase FIR Lattice Filter. 
e(k) = d(k) - y(k) 

(5.31) 

- d(k) - s M [f^lk) + gM-l(k-D-l)] ' k M l f M-l( k ' D ) + ^M-l^' 1 ^ 

where d(k) is the desired signal. A schematic of the linear phase FIR lattice realization 
is shown in Figure 5.6. From Eq.(5.31) we see that y(k) is a function of both Siyj and 
k^p But fi\j.i(k) and g^j.jlk-l) are functions of k-^pj and so on. Therefore, in order 
to minimize the mean square error, we define a cost function 



J = E { e 2 (k) } 

and a gradient element 



(5.32) 
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z(k) = 



(5.33) 



5y(k) 




Then following the treatment in Eqs.(5.6)-(5.11) the LMS algorithm for the linear 

phase lattice can be written as 



kj(k+ I) - kj(k) + 2 |i k e(k) z(k) 

s M (k+l) - s M (k) + 2 (i s e(k) (f M .,(k) + g fvl - l< k 'D - 1 >] 
where j= 1,2, From Eqs.(5.29) and (5.30), the gradient vector z(k) is given by 

z(k) - 0 MJ (k) + ■P Mij (k-D) 

- s M (k) [3> M -l,j< k ) + ^M-l.jfk D-I)] + k M (k) [® M -l.j( k -D) 

+ 'F M . I j (k-l)] + (f M . ; (k-D) - g M .,(k-l)] 6 Mij 

where 6jyj j = (5k^j/0kj), and <I> m j and 'F m j are obtained on the same lines as in 
Eqs.(5. 15)-(5.17). The resulting LMS algorithm is summarized in Table 2. 

Simulation Results : 

The performance of the algorithm summarized in Table 2 has been observed by 
computer simulation. The configuration used is the system identification experiment as 
discussed in Section (IV. D). We have used two different fixed coefficient linear phase 
FIR filter examples, one with N even and the other with N odd. They are 
H 3 (z) = 0.154 + 0.462Z* 1 + 0.462z‘ 2 + 0.154z' 3 
and 

H 4 (z) = 0.15 - 0.45z‘* + 0.36z' 2 - 0.45z' 3 + 0.1 5z'^ 

The learning curves obtained using the new algorithm are shown in Figures 5.7 to 5.10. 
Figures 5.7 and 5.S are obtained learning curves with linear phase FIR transfer 

function H 3 (z). Figures 5.9 and 5.10 are obtained learning curves with linear phase 
FIR transfer function H 4 (z). 
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TABLE 2 

LINEAR PHASE FIR LATTICE ALGORITHM 

Initialization: 

K( 0) = 0 
s jvj( 0 ) = 1 
Lattice: 

x(k) = f 0 < k ) = 

f m ( k ) = s m (k)f m . l( k) + Vlcjg^idc-l) 
g m ( k ) = s m ( k) g m _]_( k-1 ) + k m (k)f m . 1 (k) 
m = 1 , 2 , . . . , M 

y(k) = Sjyj [fj v j_^(k) + < ?iyi-i( k-D- - ) 1 + Rjyj £ ^M-l( k ~^) 

+ %l-l( k- ^^ 

Update Equations: 

k,(k+l) = kj(k) + 2 (k)] e(k) Zj(k) 

s M (k+l) = s M (k) + 2 [ ?x s /<y 2 s ( k ) ] e(k) [ f M-1 ( k) 

+ g M _i(R-D-i)] 

where ji k and n s are the adaptation constants, 

<r 2 k _.(k) = X <J 2 k _.( k-1) + (1-X) [0> M/j (k) +T M/ j(k-D)] 2 
and 

<y 2 s (k) = X cr 2 s ( k-1) + ( 1-X) [f M _ 1 (k) + g M-1 ( k-D-1) 3 2 

are estimations of power in Zj(k) and 
[ £m-i( k ) +g M-l( k ”D-l )], respectively and X is 
a positive weighting constant, 0<X^1. 

Gradient Vector Elements: 





= >P 0 /j (k) = o 


„(k) 


= s<(k)<t> H , ,(k) + k i (k)'F 1 . 1 .(k-1) 




- g i . 1 (k-l)6 i/ j 


'P i ,j(k) 


= s i (k)T i . 1/ j(k-l) + k i (k)a> i . lj (k) 




+ ^i-l^ k ^i,j 


Zj(k) = 


*M,j( k ) + T M/j (k-D) 



61 




Figure 5.7 Learning Curve (p. = 0.1). 




Figure 5.8 Learning Curve (p = 0.3). 
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Figure 5.9 Learning Curve (p. = 0.03). 



o 




Figure 5.10 Learning Curve = 0.1). 
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D. SPECTRAL ESTIMATION 



In this section, we extend the linear phase FIR adaptive lattice algorithm derived 
in the previous section to the estimation of spectral lines in white noise. The spectral 
estimation problem is some what different from the system identification experiment 
that we have been considering so far. In a typical spectral estimation problem, we are 
given a time series and we need to estimate its spectrum. A suitable configuration that 
is frequently used for this purpose is shown in Figure 5.11. Comparing this with 
Figure 4.5, we notice that we have the given time series at the input of the adaptive 
filter and the filter updates its coefficients to minimize the mean square output. In 
doing this, we assume that the input process x(k) has been generated by passing a 
white noise sequence through a filter, say I(z), and the adaptive filter attempts to 
realize an inverse filter, say H(z)=l/I(z). Considering that the adaptive filter has 
sufficient degrees of freedom, the output of the filter will be a white noise sequence. 

The adaptive algorithm summarized in Table 2 can be used in spectral estimation 
after appropriate modification. In the present case we minimize the cost function, 

J = E{y 2 (k)} 

rather than J = E{e 2 (k)}. The resulting update equations can be shown to be 



kj(k+ 1) = kj(k) - 2 [H k /<r\ (k)] e(k) Zj (k) (5.36) 

s M (k+ 1) - s M (k) - 2 [M s /<* 2 s (k)] e(k) [f M .j(k) 4- g M . 1 (k-D-l)] (5.37) 

where Zj(k), <? k .(k) and <7 s (k) are as defined in Table 2. 

Using the adaptive algorithm in Eqs.(5.36) and (5.37) we can obtain values of the 
reflection coefficients and the s^j coefficient. We then obtain the equivalent 
polynomial. B(z), by going through the standard step up procedure given by [Ref. 1]: 

°m,0 - s m 

h = k 
m,m ^m 

^m.i ~ s m ^m-l,i + ^m^m-l,m-i 
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where i= l,2,...,m-l and m= Finally the required linear phase transfer 

function is obtained as follows: 

H N .i(z) = F m (z) + z' D G m (z) 

where Fjyj(z) = B(z) and G^(z)= z'^B(z'*). The spectrum is computed as follows: 

1 

S( 0 = : : 

| h 0 + hj e'i 1(0 + ... + h N .j e-K 11 ' 1 ) 03 | 2 

where to = 27tf. and f is the normalized frequency (with respect to the sampling 
frequency) in the range, 0 :< f< 0.5. 

. Simulation Results : 

The input process x(k) consists of a signal in noise, given by 
x(k) = s(k) + w(k) 

where s(k) may be a single, or multiple sinusoids and w(k) is a zero mean unit variance 
white noise sequence. In the following we consider several examples using parameters 
ranging from a single sinusoid to 4 sinusoids, SNRs from 30dB to lOdB, and filter 
order, M, from 2 to 30. 

Example I : We consider a single sinusoid given by 



x(k) = V2 cos(2rtfk) + w(k) (5.38) 

where f=0. 15, SNR=30dB. 

Figures 5.12. 5.13. 5.14, and 5.15 are plots of frequency spectrum with different 
order of lattice M and adaptation constant ft. 

Figure 5.12 is the plot of the case M = 2 and = 0.015. We see that the peak is 
at f=0.15 and no spurious responses. 

Figure 5.13 shows the frequency spectrum of M = 4 and ft = 0.01. We observe 
that one peak is at f=0.15 and a spurious response whose magnitude is about -53dB at 
f=0.37. 
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Figure 5.11 Spectral Estimation Mode. 

Figure 5.14 shows the frequency spectrum of M= 10 and fi = 0.03. There are one 
peak at f= 0.15 and four spurious responses at f=0.05, 0.25, 0.35 and 0.45. The largest 
spurious response is about -15 dB at f=0.45. 

Figure 5.15 shows the frequency spectrum of M = 20 and n = 0.03. There are one 
peak at f=0.15 and nine spurious responses. The largest spurious response is about 
-27 dB at f= 0.125. 

Through the Figures 5.12, 5.13, 5.14, and 5.15 we can determine that the best 
model order to detect the one sinusoidal signal is VI = 2. 

Next, we fix the order of lattice M. adaptation constant u and iteration number 
k. Figures 5.16, 5.17 and 5.18 are the plots of frequency spectrum with changing the 
signal to noise ratio (SNR). 

Figure 5.16 is the plot of M = 20, )t = 0.05, k=3000 and SNR=30dB. There is a 
peak at f = 0.15 and the largest spurious response is about -27dB at f= 0.125. 
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Figure 5.17 is the plot of M = 20, n = 0.05, k=3000 and SNR=20dB. The largest 
spurious response at f= 0.425 is larger than the desired response. At this case, we 
cannot estimate the frequency of sinusoid. 

Figure 5. IS is the plot of M = 20, n = 0.05, k=3000 and SNR=10dB. Three of 
the spurious responses are larger than the desired response. The desired response is 
about -13dB. 

From Figures 5.16, 5.17 and 5.18, if the SNR is getting smaller then estimating 
sinusoidal frequency is more difficult. 

Example 2: Next, to estimate two sinusoidal frequencies in noise, set the input x(k) as 



x(k) = V2 { cos(27tfjk) + cos(27tf2k) } + w(k) (5.39) 

where the normalized frequencies of signals fj = 0.15 and f 2 = 0.25 and set SNR= 30dB. 

Figure 5.19 shows the frequency spectrum of M = 4 and u = 0.02. There are two 
peaks at fj = 0.15 and f2='0.25 and no spurious responses. 

Figure 5.20 shows the frequency spectrum of M = 20 and p = 0.022. There are 
two peaks at fj = 0.15 and f^O.25, and eight spurious responses. The largest spurious 
response is about -20dB at f= 0.375. 

From Figures 5.19 and 5.20, the best model order to estimate the frequencies of 
two sinusoidal signals is M = 4. 

Example 3: Next, to estimate four sinusoidal frequencies in noise, set the input x(k) as 



x(k) = V2 { cos(27tfjk) + cos(27tf?k) + cos(2Ttfjk) + cos(2^f^k) } + w(k) (5.40) 

where the normalized frequencies of signals fj = 0.05, f-> = 0. 15. F, = 0.25, Tj = 0.35, and 
set 3NR= 30dB. 

Figure 5.21 shows the frequency spectrum of M = 8 and |1 = 0.015. There are four 
peaks at fj = 0.07, f 2 = 0.15, fj = 0.27 and = 0.375, and no spurious responses. 

Figure 5.22 shows the frequency spectrum of M = 30 and p = 0.014. There are 
four peaks at fj=0.05, F 2 = 0.15, fj = 0.25 and f^ = 0.35, and eleven spurious responses. 

The largest spurious response is about -23dB at f=0.45. 
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From Figures 5.21 and 5.22, the best model order to estimate the frequencies of 
four sinusoidal signals is M = 8. 




Figure 5.12 Frequency Spectrum (1 sinusoid, M = 2, = 0.015). 
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Figure 5.13 Frequency Spectrum (1 sinusoid, M = 4, ji = 0.01). 




Figure 5.14 Frequency Spectrum (1 sinusoid, M= 10, n = 0.03). 
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Figure 5.15 Frequency Spectrum (1 sinusoid, M = 20, ^1 = 0.03). 




Figure 5.16 Frequency Spectrum (l sinusoid, M = 20, SNR= 30dB). 
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Figure 5.17 Frequency Spectrum (1 sinusoid, M = 20, SNR= 20dB). 




Figure 5.18 Frequency Spectrum (1 sinusoid, M = 20, SNR= 10dl3). 
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Figure 5.19 Frequency Spectrum (2 sinusoids, M = 4, p = 0.02). 




Figure 5.20 Frequency Spectrum (2 sinusoids, M = 20, |i = 0.022). 
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Figure 5.21 Frequency Spectrum (4 sinusoids, M = 8, = 0.015). 



O 




Figure 5.22 Frequency Spectrum (4 sinusoids, M = 30, n = 0.014). 
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E. SUMMARY, CONCLUSION AND SUGGESTED FUTURE WORK 

The thesis investigates the application of lattice structures in Prony's method of 
spectral line estimation. The complete solution to Prony's method consists of three 
steps: (i) representing a given process of M sinusoids in noise in terms of complex 
exponentials, (ii) finding roots of a symmetric polynomial, and (iii) estimating the 
frequency, phase and amplitude information. In this thesis, however, we have 
emphasized only the frequency estimation problem. 

The underlying theory of Prony's method has been briefly reviewed in Chapter II. 
By using a modified Prony's approach, we observed that the frequency estimation (as 
part of Prony's method) requires a polynomial with a linear phase property. The 
basics of the linear phase FIR filter and the lattice structure have been included in 
Chapter III. Also we have shown that a linear phase FIR transfer function can be 
realized from an FIR lattice using D number of additional unit delays and an adder, 
where D is determined by the order of the linear phase polynomial. 

The principles of adaptive filtering and the LMS algorithm have been briefly 
studied in Chapter IV. We have derived a new LMS algorithm for the FIR non-linear 
phase ans linear phase transfer functions in Chapter V. The problem spectral 
estimation using the new adaptive algorithm has been addressed, and the simulation 
results are included. 

The significant outcome of the work is the derivation of an LMS type algorithm 
for a linear phase lattice structure and its application to spectral line estimation as part 
of Prony's method. As has been shown, the new algorithm yielded faster convergence 
rate compared to previously reported approximate algorithms. The application of this 
algorithm to spectral estimation produced high spectral resolution as illustrated by the 
simulation results. 

However, we have observed spurious spectral responses, when the model order is 
higher that the required. Also w r e have observed that the SNR performance of the 
algorithm needs to be improved. For input SNRs less than about 10 dB, the algorithm 

performed poorly. Finally, we have dealt with only the frequency estimation part of 
Prony's method. In a future work, some of the above mentioned problems may be 
taken up. 
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APPENDIX A 

EXAMPLES OF OBTAINING A LATTICE STRUCTURE FOR THE FIR 

TRANSFER FUNCTION 



Example A. I 

Obtaining a lattice structure for the general FIR transfer function. The unit 
sample response of a FIR transfer function is given by 
H 3 (z) = 0.5 + 0.25Z* 1 + 0.125z’ 2 + 0.0625z* 3 

Solution : Given unit sample response is 

H 3 (z) = 0.5 + 0.25Z* 1 + 0.125z' 2 + 0.0625z' 3 (A.l) 



Using Eq.(3.22), we have polynomial for 3 sections 



H 3< z > = b 3,0 + b 3,l z ' 1+b 3,2 z ' 2 + b 3,3 z 



- J 



Comparing Eqs.(A.l) and (A. 2), we have 



b 3,l = °- 25 

b^ 9 = 0. 125 
b 3,3 = °- 0625 

Starting with m= 3 we have from Eq.(3.46) 



(A. 2) 



(A. 3) 



K.3 — b 3 3 ~ 0.0625 
s 3 = b 3,0 = 0t5 



(A.4) 



Now, we need to generate the coefficients for ^(z) and from Eq.(3.46) 
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c k ♦ - k h • 
k — m m.i mm.m-i 
D m-l,i s l . i. l 



s m * k m 



(A. 5) 



And for m = 3 and i = 0 we have 

h _ s 3 b 3.0 * k 3 b 3.3 

° 2,0 Tr~T^ 



b 2,0 “ 



(0.5X0.5) - (0.0625)(0.0625) 
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W 



(0.0625)" 



= 1 



(A. 6) 



Next, for m = 3 and i = 1 we get 



5 2,1 



= s 3 b 3.1 ' k 3 b 3.2 



so - ki- 



(0.5X0.25) -(0.0625X0.125) __ 

bo i = , j = 0.47619 



J 2,l 



(0.5) - (0.0625)" 



(A.7) 



and finally 



b 2,2 ” 



s 3 b 3.2 * k 3 b 3.1 



So - ki 



^ _ (0.5)(0.125) - (0.0625)(0.25) _ a 

b 7 o — t j — 0.19048 

l ' L (0.5) 2 - (0.0625) 2 



from Eq.(3.46), we have 



k 2 “ b 2.2 = 0.19048 

s 2 = b 2,0 = 1 

The new polynomial is 

H 2 (z) = b 2>0 + b 21 z* 1 + b 2)2 z* 2 



(A.8) 



(A. 9) 
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(A. 10) 



H 2 (z) = 1 + 0.47619Z' 1 + 0.19048Z' 2 



for m= 2 and i= 0, we have 



b 



1.1 “ 



b 2.0 * k 2 b 2.2 
1 - 



1 -(0.19048X0.19048) 

1 - (0.19048) 2 



(A.l 1) 



and finally 



b 



1,1 " 



b 2.1 * k 2 b 2.1 
l * V 



0.47619 - (0.19048)(0.47619) 

b, , = --t — = 0.40000 

kl 1 - (0.19048)" 



From Eq.(3.46), we have 



(A. 12) 



k l = b l,l = 0- 40000 
S 1 = b l,0 = 1 



(A. 13) 



Therefore, reflection coefficients and s coefficients of the FIR lattice structure are 



k 

k 

k 

s 



1 

2 

3 

1 



So — 



s 3 = 



0.40000 

0.19048 

0.0625 

1 

1 

0.5 
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Example A. 2 (N= 4:everi) 

Obtaining a lattice structure for the linear phase FIR transfer function. A linear 
phase FIR transfer function is 

H 3 (z) = 0.154 - 1 - 0.462Z’ 1 + 0.462z‘ 2 + 0.1 54z* 3 (A. 14) 

Solution : Eq.(A.14) can be written in the form of Eq.(3.35) 

] 0 1 

H 3 (z) = aQ + ajZ‘ A -f z“ (aj + aQZ _i ) (A.15) 

Comparing Eqs.(A.14) and (A. 15), we get the following relationships. 



aQ = 0.154 
aj = 0.462 



(A. 16) 



To determine the lattice reflection coefficients of the corresponding linear phase FIR 
lattice structure and the number of unit delay, we use Eq.(3.29) which is 

H N .,(z) = F m (z) + z- D G m (z) (A. 17) 



From Eq.(A.14), we have N = 4(even), the order of the polynomials Fj^(z) and G-yj(z), 
M = (N/2)-l = (4/2)-l = 1, and the number of unit delays, D = N/2 = 2. Now, 
from Eq.(3.37), we may have the forward prediction error transfer function as follows: 



F l(z) = a 0 + aj z' 1 

= 0.154 + 0.462 z* 1 



(A.18) 



for M= 1, Eq.(A.i8) can be rewritten as 

F,(z) = b, 0 + bj^z' 1 



(A. 19) 
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From Eqs.(A. 18) and (A. 19), we get the desired reflection coefficient and s coefficient 
of the linear phase FIR lattice structure as follows: 

k l = b l,l = °* 462 
S 1 = b l,0 = °- 154 

Example A. 3 (N = 5:odd) 

Obtaining a lattice structure for the linear ^phase FIR transfer function. A linear 
phase FIR transfer function is 

H 4 (z) = 0.15 - 0.45Z' 1 + 0.36z‘ 2 - 0.45z‘ 3 + 0.15z‘ 4 (A.20) 

Solution : Eq.(A.20) can be written in the form of Eq.(3.30) 

H 4 (z) = aQ + ajZ** + (l/2)a2Z*^ + z' 2 {(l/2)a2 + ajz** + aQZ* 2 } (A.21) 
Comparing Eqs.(A.20) and (A.21), we get the following relationships. 



To determine the lattice reflection coefficients of the corresponding linear phase FIR 
lattice structure and the number of unit delay, we use Eq.(3.29) which is 



From Eq.(A.20), we have N=5(odd), the order of the polynomials Fj^(z) and G-yj(z), 
M = (N-l)/2 = (5-l)/2 = 2. and the number of unit delays. D = (N-l)/2 = 2. Now. 

from Eq.(3.32), we may have the forward prediction error transfer function as follows: 



aQ = 0.15 
aj = -0.45 
a2 = 0.36 



(A. 22) 



h n-i( z ) = f m( z ) + Z * D g m( z ) 



(A.23) 



F 2 (z) = aQ + aj z*^ + (1/2) a 2 z* 2 



(A. 24) 



0.15 - 0.45 z* 1 + 0.18 z* 2 



for M = 2, Eq.(A.24) can be rewritten as 
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(A. 25) 



F 2< z > = b 2,0 + b 2,l z '‘ + b 2,2 z ' 2 

From Eqs.(A.24) and (A.25), we have 
k .2 ~ b2 2 = 0.18 
s 2 = b 2,0 = °‘ 15 

Now, we need to generate the coefficients for Fj(z) and from Eq.(3.46) 



b l,0 " 



S2b 



2,0 

T ~ 



- k->b 



12 . 



- k- 



(0.15>(0.15) - (0.18)(0.18) 
0.15 2 - 0.18 z 



(A. 26) 



b l,l = 



s 2 b 2.1 I k 2 b 2.1 



- k- 



(0.15)(-0.45) - (0.18)(-0.45)_ 
0.15 2 - 0.18" 



-1.36364 



(A. 27) 



From Eqs.(A.26) and (A.27), we have 
ki = bi i = -1.36364 
S 1 = b l,0 = 1 

Therefore, reflection coefficients and s coefficients of the linear phase FIR lattice 
structure are 

= -1.36364 
= 0.18 
= 1 
= 0.15 



S 1 

s 2 
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APPENDIX B 

COMPUTER PROGRAMS 



rogram l* ******************* ********** xx ************ 7 ' : ** *********** 

* This program is designed for the system identification experiment. 

* which is shown in Section (IV. D). The learning curves can be obtained* 

* by plotting the error, E, versus the update iteration, k. * 



* 

* 

* 

k 

X 

k 

* 

* 

X 

* 

* 

A 

* 

* 

"k 

k 

* 

* 

* 

■k 

k 

k 



k 



k 



k 



Variable Definition 

ISEED : seed for the random number generation (white noise) 
AMU : adaptation constant 

k : time index 

M : order of the FIR transfer function or total number of 

lattice sections 
NB s number of iterations 

A(I) : coefficients of the FIR transfer function 

B (I ) : reflection coefficients of the lattice 

F(I) : forward prediction error 

G(I) : backward prediction error 

GD(I) : delayed backward prediction error 

SGL(I): estimations of power 

W(K) : input of both FIR and lattice 

YF(K) : output of the FIR filter 

YL(K) : output of the lattice filter 

ER(K) : squared error 

Variable Declaration 

INTEGER ISEED, K,I,J,N,M,NB 

REAL A(100) ,B(100) ,F(100) ,G(100) ,GD(100) ,YF(10000) ,YL(10000) 
REAL X(100) ,SGL(100) ,ER(10000) ,W(10000) , Y(10000) ,AMU,E 

Set Adaptation Constant U and Number of Iterations NB 



WRITE (5,1) 

1 FORMAT ( 5X , 1 AMU 1 , 5X , 1 NB 1 ) 

READ (5 , *) AMU,NB 

Initialization 



DO 10 K=1 , 100 
A(K)=0 
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B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
SGL(k)=l .0 
10 continue 

DO 15 K=l, 10000 
YF (K)=0 
YL(K)=0 
ER(K)=0 
W(K)=0 
Y(k)=0 

15 CONTINUE 
E=0 . 

R=1 

ISEED=343169 

M=2 

A(l) = 1.0 

A(2) = -0.89 
A(3) = +0.25 

* 

* Random Number Generation 

(mean zero, unit variance, white sequence) 

* 

DO 20 N=1,NB 

CALL LNORM(ISEED,RN, 1,1,0) 

W(N) = RN 
20 CONTINUE 

* 

* FIR Filter Calculation 

* 

DO 30 K = 1 ,NB 
X(1)=W(K) 

YF (K)= A(1)*X(1) 

DO 40 1=1, M 

YF(K)=YF(K)+A(I+1)*X(I+1) 

40 CONTINUE 

DO 45 1=1, M 

X(M+2-I)=X(M+l-I) 

45 CONTINUE 

30 CONTINUE 

* 

* Lattice Filter Calculation 

* 



DO 50 K =1 ,NB 
F(1)=W(K) 
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G(1)=W(K) 

DO 60 I = 1 ,M 

F(l+1)=F(I)+B(I)*GD(I) 

G(I+1 )=GD (I )+B (I )*F( I ) 

60 CONTINUE 

YL(K)=F(M+1) 

* 

* Calculating the Error 

* 

E = YF(K) - YL(K) 

k 

* Updating the Reflection Coefficients 

* 

CALL LMS ( B , GD , E , AMU , SGL , M ) 

k 

ER(K)=E**2 
DO 70 J=1 ,M 
GD( J)=G( J) 

70 CONTINUE 

IF (K.EQ.R) THEN 

WRITE (6 , 300) K, 3(1) ,3(2) ,ER(K) 

R=R+50 
END IF 
Y(k)=E 

50 CONTINUE 

300 FORMAT (3X, 16 , 5X, 3 (F10 . 7 ,4X) ) 

* 

* Plotting the Learning Curve 

* 

CALL PLOT(Y,N) 

STOP 

END 

* 

* 

k 

SUBROUTINE LMS ( B , GD , E , AMU , SGL , M ) 

REAL B(100) , GD(100) , SGL (100) ,E,AMU 
INTEGER M 

DO 200 1=1 ,M 

SGL( I)=0 . 9*5GL(1 )+0 . 1*(GD(I) *GD (I) ) 

200 CONTINUE 

DO 210 1=1 ,M 

B(I)=B(I)+(AMU/SGL(I))*E*GD(I) 

210 CONTINUE 

RETURN 
END 

* 
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k 

k 



10 

c 

c 



c 



SUBROUTINE PLOT(Y,N) 

DIMENSION Y(N) ,X(10000) 

ISTP=N/10 
DO 10 J=1,N 
X(J)=J 
CALL TEK618 
CALL PRTPLT(72 ; 6) 

CALL SHERPA( 'PARKPARK' , 'A' ,3) 

CALL PHYSOR(l . , 1 . ) 

CALL RESET ('ALL') 

CALL PAGE (8. 5 ,11.0) 

CALL HWROT ( ' AUTO 1 ) 

CALL XINTAX 

CALL AREA2D( 5 .0,2.8) 

CALL HEIGHT (0.10) 

CALL COMPLX 

CALL SHDCHR(90. 0,1, 0.002,1) 

CALL HEADIN( 'LEARNING CURVE$ ' , 100 , 2 . 0 , 1 ) 

CALL XNAME ( ' ITERATIONSS ' , 100) 

CALL YNAME ( ' ERROR$ 1 ,100) 

CALL MESSAG( 'ADAPTIVE LATTICE (AMU= . 5 ,FIG4 . 8) $ 1 , 100 , 3 . 0 , -0 . 8) 
CALL FRAME 

CALL GRAF(0 , ISTP ,N, -3. 0,1. 5, 3.0) 

CALL THKCRV(0 .02) 

CALL CURVE (X,Y,N,0) 

CALL ENDPL(O) 

CALL DONEPL 

RETURN 

END 
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* This program is designed for the system identification experiment * 

* using the LMS algorithm which was derived in Section (V.B). * 

* 

* 

* 

INTEGER ISEED,K, I , J,N-,M,NB 

REAL A(100) ,B(100) ,F(100) ,G(100) ,GD(100) ,x(100) ,YF(5000) , YL(5000) 

REAL ER(5000) ,W(5000) ,SGL(100) ,PH(100 , 100) ,PS (100 , 100) 

REAL PSD(100 , 100) ,GR(100) ,Y(5000) ,AMU,E 

3*C 

* Set Adaptation Constant p. and Number of Iterations NB 

* 

WRITE (5,1) 

1 FORMAT ( 5K , 1 AMU 1 , 5X , 1 NB ' ) 

READ (5,*) AMU , NB 

* 

* Initialization 

* 

E=0 . 

R=1 

DO 5 K= 1,5000 
W(K)=0 
YF(K)=0 
YL(K)=0 
ER(K)=0 
Y (K)=0 

5 CONTINUE 

DO 10 K=1 ,100 
A(K)=0 
B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD (K)=0 
SGL(K)=1 .0 
GR( K)=0 

10 CONTINUE 

DO 15 K=1 , 100 

DO 15 L=1 , 100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 

15 CONTINUE 

ISEED=343169 

M=2 
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A(l) = 1.0 

A(2) = -0.89 
A(3) = +0.25 

Random Number Generation 
(mean zero, unit variance, white sequence) 

DO 20 'N=l ,NB 

CALL LNORM(ISEED,RN, 1,1,0) 

W(N) = RN 
20 CONTINUE 

FIR Filter 

DO 30 K = 1 ,NB 
X(1)=W(K) 

YF (K)= A(1)*X(1) 

DO 40 1=1, M 

YF(K)=YF(K)+A(I+1)*X(I+1) 

40 CONTINUE 

DO 45 1=1, M 

X(M+2-I)=X(M+l-I) 

45 CONTINUE 

30 CONTINUE 

Lattice Filter 

DO 50 K =1 ,NB 
F(1)=W(K) 

G(1 )=W(K) 

DO 60 I = 1 ,M 

F(I+1)=F(I)+B(I) *GD ( I ) 

G(I+1 )=GD( I )+B (I )*F ( I ) 

60 CONTINUE 

YL(K)=F(M+1) 

E = YF(K) - YL(K) 

Updating the Reflection Coefficients 

CALL MLMS (?H,?S,?SD,GR,F,G,GD,3,SGL,Z,AMU,M) 

ER(K)=E**2 
DO 70 J=1,M 
GD(J)=G(J) 

70 CONTINUE 

IF (K.EQ.R) THEN 

WRITE (6 , 300 ) K, B(l) ,B(2) ,B(3) ,B(4) ,B(5) ,B(6) ,B(7) ,B(8) ,E 
R=R+30 
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END IF 
Y (K)=E 

c WRITE (6 , 100 ) K,W(K) ,YF(K) ,YL(K) ,E,ER(K) 

50 CONTINUE 

300 FORMAT ( IX, I6,2X,9(F10.7,2X)) 

100 FORMAT (3X , 13 , 2X , 5 (F10 . 7 , 2X) ) 

* 

* Plotting the Learning Curve 

* 

CALL PLOT(Y,N) 

STOP 

END 

A 

* 

* 

SUBROUTINE MLMS (PH, PS ,PSD ,GR, F ,B ,BD , R, SGL,EK, AMU,N) 

DIMENSION PH( 100 , 100 ) ,PS(100,100) , PSD (100 , 100 ) , F( 100 ) , B ( 100) 
1 ,BD(100) ,R(100) , GR(100) , SGL(lOl) 

GR(N)=BD(N) 

DO 200 1=2, N 
?H(I,1)=BD(N+1-I) 

PS(I , 1 )=F (N+l-I) 

IT=I-1 

DO 10 K=1 , IT 

PH (I ,K+1)=PH(I , K)+R(N+1-I+K)*PSD(I ,K) 

10 PS (I ,K+1 )=PSD (I , K)+R(N+1-I+K) *PH(I , K) 

DO 20 K=1 , IT 
20 PSD(I,K)=PS(I,K) 

200 CONTINUE 

DO 210 K=2,N 

210 GR(N+1-K)=PH(K,K) 

DO 211 K=1,N 

211 SGL(K)= . 90*SGL(K)+ . 10*GR(K)*GR(K)+1 . 0 
DO 220 1=1, N 

R(I)=R(I)+(AMU/SGL(I))*EK*GR(I) 

IF(R(I) .GE.1.0) R(I )=0 . 

IF(R(I).LE.-1.0) R(I )=0 . 

220 CONTINUE 
RETURN 
END 

* 

* 

* 

SUBROUTINE PLOT(Y,N) 

DIMENSION Y(N) ,X(5000) 

ISTP=N/10 
DO 10 J=1,N 
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10 

c 

c 



c 



X(J)=J 

CALL TEK618 

CALL PRTPLT(72,6) 

CALL SHERPA ( 1 PARKPARK ' , ' A ' , 3 ) 

CALL PHYSOR(l . , 1 . ) 

CALL RESET ('ALL 1 ) 

CALL PAGE (8.5,11.0) 

CALL HWROT ( 1 AUTO 1 ) 

CALL XINTAX 

CALL AREA2D(5 .0,2.8) 

CALL HEIGHT(O.IO) 

CALL COMPLX 

CALL SHDCHR(90. 0,1, 0.002,1) 

CALL HEADIN( 'LEARNING CURVE$ ',100 , 2 . 0 , 1 ) 

CALL XNAME ( 1 ITERATIONS$ 1 , 100 ) 

CALL YNAME( 1 ERROR$ ' ,100) 

CALL MESSAG( 'ADAPTIVE LATTICE (AMU= . 5 , FIG5 . 5) $ 1 , 100 , 3 . 0 , -0 . 8) 
CALL FRAME 

CALL GRAF (0 , ISTP ,N, -2 .5,1.25,2.5) 

CALL THKCRV(0.02) 

CALL CURVE (X,Y,N,0) 

CALL ENDPL(O) 

CALL DONEPL 

RETURN 

END 
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kkkkkkp fo^r^rn 3 

* This program is designed for the system identification experiment. * 

* The LMS algorithm shown in Section (V.C) was extended to the linear * 

* phase FIR lattice filter. * 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

k 

•k 

k 

INTEGER ISEED,K,I, J ,N,M,NB , R ,MA,ML 

REAL A(100) ,B(100) ,F(100) ,G(100) ,GD(100) ,YF(3000) ,YL(3000) 

REAL X(100) ,H(100) ,YD(3000) ,ER(3000) ,W(3000) ,YFF(300Q) ,YLL(3000) 

REAL SGL(IOO) ,PH( 100 , 100 ) , PS (100 , 100) ,PSD(100 , 100) 

REAL GR(100) ,GRD( 100 , 100 ) ,GRB( 100) ,Y(3000) ,AMU,E,C 

* 

* Set Adaptation Constant and Number of Iterations NB 

* 

WRITE (5,1) 

1 FORMAT (5X, 'AMU 1 , 5X, 'NB' ) 

READ (5,*) AMU , NB 

* 

* Initialization 

k 

E=0 . 

R=1 

M=4 

ISEED=343169 
DO 10 K=1 , 100 
A(K)=0 
B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
H(K)=0 
SGL(K)=1 . 0 
GR(K)=0 
GR3 (K)=0 
10 CONTINUE 

DO 11 X= 1,3000 
W ( K ) =0 
YF(K)=0 
YL(K)=0 
YFF(K)=0 
YLL(K)=0 
YD(K)=0 
ER(K)=0 
Y (K)=0 
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11 CONTINUE 

DO 15 K=1 , 100 
DO 15 L=1 , 100 
?H(K,L)=0 
PS (K,L)=0 
PSD(K,L)=0 
GRD(K,L)=0 
15 CONTINUE 

c C= . 154 

c A(l)=. 154/C 

c A(2)=. 462/C 

C= . 15 

A(1 )= . 15/C 
A(2)=- .45/C 
A(3)=. 36/ C 

A 

* Random Number Generation 

* (mean zero, unit variance, white sequence) 

A 

DO 20 N=1 ,NB 

CALL LNORM(ISEED,RN, 1,1,0) 

W(N) = RN 

20 CONTINUE 

IF (MOD(M,2).EQ.O) MA=M/2 
IF (MOD(M,2).NE.O) HA=(M+l)/2 
IF (MOD(M,2).EQ.O) ML=MA 
IF (MOD(M,2).NE.O) ML=MA-1 

A 

* Separation of Coefficients 

A 

IF (M0D(M,2) .EQ.0) THEN 
DO 21 1=1, MA 
H ( I ) =A ( I ) 

H(MA+2+I)=A(NA+l-I) 

21 CONTINUE 

H(MA+1 )=A(MA+l)/2 
H (MA+2 ) =A (MA+1 ) / 2 
END IF 

IF (MOD (M , 2 ) . NE . 0 ) THEN 
DO 22 1=1, MA 
H(I)=A(I) 

H(HA+I)=A(MA+1-I) 

22 CONTINUE 
END IF 

A 

* LINEAR PHASE FIR FILTER 
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DO 



40 



41 



42 



43 



45 

30 



* 

* 



60 



* 

■k 

■k 
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30 K = 1 ,NB 
X(1)=W(K) 

YF(K)=H(1)*X(1) 

IF (MOD(M,2).EQ.O) THEN 
DO 40 1=1, MA 

YF(K)=YF(K)+H(I+1)*X(I+1) 
CONTINUE 

YFF (K)=YF(K) 

DO 41 1=1 ,MA+1 

YFF(K)=YFF(K)+H(MA+1+I)*X(MA+I) 
CONTINUE 
END IF 

IF (MOD(M,2) .NE.O) THEN 
DO 42 1=1 ,MA-1 

YF(K)=YF(K)+H(I+1)*X(I+1) 
CONTINUE 

YFF (K)=YF (K) 

DO 43 1=1, MA 

YFF(K)=YFF(K)+H(MA+I)*X(MA+I) 
CONTINUE 
END IF 

DO 45 1=1, M 

X(M+2-I)=X(M+l-l) 

CONTINUE 
CONTINUE 

LATTICE FILTER 



DO 50 K =1 ,NB 
F(1)=W(K) 

G(1)=W(K) 

DO 60 I = 1,ML 

F(I+1)=F(I)+B(I)*GD(I) 
G(I+1)=GD(I)+B(I)*F(I) 
CONTINUE 
YL(K)=F (ML+1) 

YD (K)=G(ML+1 ) 
YLL(K)=YL(K)+YD(K-MA) 

Z = YFF (K) - YLL(K) 



Updating the Reflection Coefficients 

CALL MLMS ( P , Q , QD , GRB , GRD , F , G , GD , B , SGL , E , AMU , ML , MA ) 
ER(K)=E**2 
DO 70 J=1,ML 
GD( J)=G(J) 

CONTINUE 
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Y(K)=E 

WRITE(6,300) K, Y(K) 

50 CONTINUE 

* 

* Plotting the Learning Curve 

* 

CALL PLOT(Y,N) 

300 FORMAT (3X, 16 , 3X,F10 . 5) 

STOP 

END 

* 

* 

* 

SUBROUTINE MLMS ( PH , PS , PSD , GRB , GRD , F , B , BD , R , SGL , EK , AMU , N , MA ) 
DIMENSION PH(100 ,100) ,PS(100,100) ,PSD(100 , 100) ,GRD(100 , 100) 
1, GRB (100) ,F(100) ,B(100) ,BD(100) ,R(100) ,GR(100) ,SGL(101) 
GR(N)=BD (N) 

GRB(N)=F (N) 

DO 200 I=2,N 
PH (I ,1)=BD(N+1-I) 

PS (I , 1 )=F (N+l-I ) 

IT=I-1 

DO 110 K=1 , IT 

PH (I ,K+1)=PH(I ,K)+R(N+1-I+K)*PSD(I , K) 

110 PS (I , K+l )=PSD( I ,K)+R(N+1-I+K)*PH(I ,K) 

DO 120 K=1 , IT 
120 PSD (I , K)=PS (I ,K) 

200 CONTINUE 

DO 210 K=2,N 
GR(N+1-K)=PH(K , K) 

210 GRB(N+1-K)=PS(K,K) 

DO 220 K=1,N 
DO 220 L=1 ,MA 

220 GRD(K,MA+2-L)=GRD(K,MA+l-L) 

DO 230 K=1,N 
230 GRD(K, 1)=GRB(K) 

DO 240 K=1 ,N 

240 GRB(K)=GRD(K,MA+1 ) 

DO 250 K= 1 , N 
250 GR(K) =GR(K) +GR3 (K) 

DO 260 K=1 , N 

260 SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 

DO 270 1=1, N 

R(I)=R(I)+(AMU/SGL(I))*EK*GR(I) 
c IF(R(I) .GE.1.0) R(I)=0. 

c IF(R(I) .LE.-1.0) R(I)=0 . 

270 CONTINUE 
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* 

* 

* 



10 

c 

c 



c 



RETURN 

END 



SUBROUTINE PLOT(Y,N) 

DIMENSION Y(N) ,X(3000) 

ISTP=N/10 
DO 10 J=1,N 
X(J)=J 
CALL TEK618 
CALL PRTPLT(72 ,6) 

CALL SHERPA ( 1 PARKPARK ' , ' A 1 ,3) 

CALL PHYS0R( 1 . , 1 . ) 

CALL RESET ('ALL') 

CALL PAGE (8.5,11.) 

CALL HWROT ( 1 AUTO 1 ) 

CALL XINTAX 

CALL AREA2D( 5 .0,2.8) 

CALL HEIGHT (0.10) 

CALL COMPLX 

CALL SHDCHR(90. 0,1, 0.002,1) 

CALL HEAD IN ( 'LEARNING CURVES ', 100 , 2 . 0 , 1 ) 

CALL XNAME ( 1 ITERATIONS $ 1 ,100) 

CALL YNAME ( 1 ERRORS 1 , 100 ) 

CALL MESSAG( 'ADAPTIVE G-M LATTICE (F5 . 10 , AMU=0 .1)$' ,100,3.0,-0.8) 
CALL FRAME 

CALL GRAF(0 , ISTP ,N , -8. 0,4. 0,8.0) 

CALL THKCRV(0 . 02) 

CALL CURVE (X,Y,N,0) 

CALL ENDPL(O) 

CALL DONEPL 

RETURN 

END 
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rogram 4 ******************************** ******************* : * ****** 

* This program is designed for the estimation of spectral lines in * 

* white noise. The input process x(k) consists of a signal in noise, 

* and a signal may be a single, or multiple sinusoids. The algorithm * 

* was derived in Section (V.D). * 

x 

* 

* 

INTEGER ISEED,K,I , J,N,M,NB,R,MA,ML,D 

REAL A(100) ,B( 100) ,F(100) ,G(100) ,GD(100) ,YL(5000) 

REAL Y(100) , YD(5000) ,W(5000) , YLL(5000) ,YD(5000) , INP(5000) 

REAL SGL( 101 ) ,PH( 100 , 100 ) , PS (100 , 100) , PSD( 100 , 100) 

REAL GRD (100 , 100) ,GRB(100) ,GR(100) 

REAL RE(100) ,IM(100) ,AJ(100) ,SD,SNR,AVG,AMP,AMU,E 

* 

* Set Adaptation Constant ^ and Number of Iterations NB 

* 

WRITE (5,1) 

FORMAT (5X, 'AMU' ,5X, 'NB' ) 

READ ( 5 , *) AMU , NB 

* 

* Initialization 

X 

c ISPEC=1000 

SNR=30 . 

SD=10**(- (SNR/20)) 

AMP=SQRT(2. ) 

AVG=0 . 

E=0 . 
c R=1 

PI=3. 141592654 
M=8 

MP1=M+1 

ISEED=343169 

IF (MOD(M,2) .EQ.0) THEN 

MA=M/2 

ML=MA 

ELSE 

HA=(M+i)/2 

ML=MA-1 
END IF 

DO 10 K=1 ,100 
A(K)=0 
B(K)=0 
F(K)=0 
G(K)=0 
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GD(K)=0 
SGL(K)=1 . 0 
GR(K)=0 
GRB (K)=0 
RE (K)=0 
IM(K)=0 
AJ(K)=0 

y(k)=o 

10 CONTINUE 

DO 11 K=1 , 5000 
W(K)=0 
YL(K)=0 
YLL (K)=0 
YD(K)=0 
ER(K)=0 
INP (K)=0 

11 CONTINUE 

DO 15 K=1 , 100 
DO 15 L=1 , 100 
PH(K,L)=0 
PS (K,L)=0 
PSD(K,L)=0 
GRD ( K , L ) =0 

15 CONTINUE 

"k 

* Random Number Generation 

* (mean zero, unit variance, white sequence) 

k 

DO 20 N=1,NB 

CALL LNORM(ISEED,RN, 1,1,0) 

W(N) = SD*RN+AVG 

20 CONTINUE 

k 

* LATTICE FILTER 

k 

DO 50 K =1 ,NB 
AK=K-1 

c INP(K)=AMP*C0S(2*?I* . 15*AK)+W(K) 

c INP(K)=AMP*<COS(2*?I*.I5^AK)+COS(2 7 ' : PI*.25*AK) )+W(K) 

INP ( K ) =AMP* ( COS ( 2*PI * . 0 5 *AI< ) +C0S ( 2*PI * . 1 5 *AK ) +C0S ( 2*?I * . 25*AK 

*)+C0S (2*PI* . 35*AK) )+W(K) 

F(1)=INP(K) 

G(1)=INP(K) 

DO 60 I = 1,ML 

F (1+1 )=F(I )+B (I )*GD(I ) 

G(I+1)=GD(I)+B(I)*F(I) 

60 CONTINUE 
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YL(K)=F (ML+1 ) 

YD(K)=G(ML+1) 

YLL(K)=YL(K)+YD(K-MA) 

E = YLL(K) 

Updating the Reflection Coefficients 

CALL MLMS (PH, PS ,PSD ,GR ,GRB ,GRD , F ,G,GD , B , SGL ,E , AMU,ML ,MA) 
ER(K)=E**2 
DO 70 J=1 ,HL 
GD(J)=G(J) 

70 CONTINUE 

c IF (K.NE.ISPEC) GO TO 50 

IF (K.EQ.NB) THEN 
c WRITE (6 , 300 ) K, INP(K) ,E 

c WRITE(6 ,301) K,B(1) ,B(2) ,B(3) ,B(4) ,B(5) ,B(6) ,B(7) ,B(8) ,B(9) 

c R=R+100 

Determine the FIR Coefficients from the Lattice Reflection 
Coefficients 

CALL STEPUP (A,B,ML) 

IF (M0D(M,2) .EQ.0) THEN 
DO 80 1=1, M/2 

80 A(M+2-I)=A(I) 

A(M/2+l )=2*A(M/ 2+1 ) 

ELSE 

DO 81 1=1 , (M+l)/2 

81 A(M+2-I)=A(I) 

END IF 

c WRITE (6,600) (I ,A(I) , 1=1 ,MP1) 

Calculate the Power Spectrum 

CALL SPEC (RE , IM, A , AJ ,M, PI ) 

D=100 

c WRITE (6,601) (.005*Z,AJ(Z),Z=1,D) 

Plotting the Spectrum 



CALL PLOT(AJ,D) 
c ISPEC=ISPEC+1000 
END IF 

50 CONTINUE 

c WRITE (6,300) (YLL(K) ,K=1 ,NB) 

300 FORMAT ( 10 (FI 0 . 7 , IX) ) 

301 FORMAT ( IX, I5,1X,10(F10.7,1X)) 
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600 FORMAT (IX, 15 , 5X, F15 .7 ) 

601 FORMAT (IX, F5 . 3 , 5X, F15 . 7 ) 
STOP 

END 

* 

* 



SUBROUTINE MLMS (PH, PS ,PSD ,GR,GRB ,GRD ,F ,B ,BD ,R, SGL ,EK, AMU,N,MA) 
DIMENSION PH ( 100 , 100 ), PS (100 , 100) ,PSD(100 , 100) ,GR(100) ,GRB(100) 
1 ,GRD(100 ,100) , F (100 ) ,B(100) ,BD(100) ,R(100) ,SGL(101) 

GR(N)=BD(N) 

GRB(N)=F (N) 

DO 200 1=2, N 
PH (I ,1)=BD(N+1-I) 

PS (I , 1 )=F (N+l-I ) 

IT=I-1 

DO 110 K=1 , IT 

PH (I ,K+1)=PH(I ,K)+R(N+1-I+K)*?SD(I ,K) 

PS (I ,K+1)=PSD(I , K)+R(N+1-I+K)*PH(I ,K) 

110 CONTINUE 

DO 120 K=1 , IT 
120 PSD (I ,K)=PS(I,K) 

200 CONTINUE 

DO 210 K=2,N 
GR(N+1-K)=PH(K,K) 

GRB(N+1-K)=PS(K,K) 

210 CONTINUE 

DO 220 K=1,N 
DO 220 L=1,MA 

220 GRD (K,MA+2-L)=GRD (K,MA+1-L) 

DO 230 K=1,N 
230 GRD(X,1)=GRB(K) 

DO 240 K=1,N 
240 GRB(K)=GRD(K,MA+1) 

DO 250 K=1 ,N 
250 GR(K)=GR(K)+GRB(K) 

DO 260 K=1,N 

260 SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 

DO 270 1=1, M 

R(I)=R(I)-(AMU/SGL(I))*EK*GR(I) 

IF(R(I) .GE.1.0) R(I)=0. 

IF(R(I).LE.-1.0) R(I )=0 . 

270 CONTINUE 
RETURN 
END 



* 

* 
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SUBROUTINE STEPUP (A,B ,ML) 

DIMENSION A(100) ,C(100) ,B(100) 

A(l)=l. 

A(2)=B(1) 

DO 30 MINC=2 ,ML 
DO 10 J=1 ,MINC 
JB=MINC- J+l 
10 C(J)=A(JB) 

DO 20 IP=2 ,MINC 

20 A(IP)=A(IP)+B(MINC)*C(IP-1) 

A(MINC+1)=B(MINC) 

30 CONTINUE 
RETURN 
END 

k 

k 

k 

SUBROUTINE SPEC(RE , IM, A, AJ ,M, PI ) 

REAL RE(100) , IM(100) ,A(100) ,AJ(100) ,Y(100) 

RE ( 1 ) =A ( 1 ) 

IM(1)=0 . 

DO 91 J=l,100 
DO 92 1=1, M 

RE(I+1)=RE(I)+A(I+1)*COS(2*I*PI*.5*J/100) 

IM(I+1)=IM(I)+A(I+1)*SIN(2*I*PI*.5*J/100) 

92 CONTINUE 

AJ ( J )=- 10 . *ALOG10 (RE (M+l ) *RE (M+l ) +IM (M+l ) *IM(M+1 ) ) 
91 CONTINUE 
TEMP=A J ( 1 ) 

DO 93 L=2 , 100 

IF (AJ(L) .GT.TEMP) TEMP=AJ (L) 

93 CONTINUE 

DO 94 L=1 , 100 
AJ(L)=AJ(L) -TEMP 

94 CONTINUE 
RETURN 
END 



SUBROUTINE PLOT(Y,N) 
DIMENSION Y(N) ,X( 100) 
c ISTP=N/10 

c DO 10 J=1,N 

CIO X(J)=J 
DF=. 5/N 
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10 



20 

c 

c 



c 



c 



K(1)=0 
DO 10 K=2,N 
X(K)=X(K-1)+DF 
XHIN=X(1) 

XMAX=X(N) 

XSTP=10*DF 
IYMIN=Y (1 ) 

IYMAX=Y(1) 

DO 20 K=2 ,N 

IF(Y(K).GT. iymax) IYMAX=Y(K) 

IF(Y(K) .LT.IYMIN) IYMIN=Y(K) 

CONTINUE 

IYSTP=(IYMAX-IYMIN)/5 

CALL TEK618 

CALL PRTPLT(72,6) 

CALL SHERPA ( 1 PARKPARK 1 , 1 A ' , 3 ) 

CALL RESET ('ALL') 

CALL PAGE (8. 5, 11.0) 

CALL HWROT ( 1 AUTO 1 ) 

CALL XINTAX 

CALL AREA2D(5 .0,2.3) 

CALL HEIGHT (0.10) 

CALL COMPLX 

CALL SHDCHR(90.0, 1,0.002,1) 

CALL HEADIN( 'FREQUENCY SPECTRUM$ 1 , 100 , 2 . 0 , 1 ) 

CALL XNAME ( 1 FREQUENCY $ 1 , 100) 

CALL YNAME( 1 MAGNITUDE (DB) $ 1 ,100) 

CALL MESSAG( 'FIGURE 5.21 (M=8 , SNR=30DB) $ ' , 100 , 3 . 0 , -0 . 8) 

CALL MESSAG( 'MODEL ORDER SELECTION^ SINUSOIDS )$', 100 , 3 . 0 , -0 . 8) 
CALL THKFRM(0 .03) 

CALL FRAME 

CALL GRAF (XMIN , XSTP , XMAX , IYMIN , IYSTP , IYMAX) 

CALL THKCRV(0.02) 

CALL CURVE (X,Y,N,0) 

CALL ENDPL(O) 

CALL DONE PL 

RETURN 

END 
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